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ABSTRACT 


The problem of a laminar water jet impinging on an 
ice surface is studied analytically. The jet flow is 
divided into two parts, a potential flow solution and a 
boundary layer solution. The potential flow problem is 
formulated as a variational statement and is solved by the 
finite element method. From the potential flow solution, 
the velocity distribution on the impingement surface is 
used as the free stream velocity in the boundary layer 
analysis. The boundary layer is solved by a Karman-Pohl- 
hausen integral technique. 

To verify the solution method the results for a two- 
dimensional jet impinging on a flat surface without melting 
are compared with exact solutions for the potential flow 
and a numerical calculation of the boundary layer. 

Solution of the problem of an axisymmetric jet imping- 
ing on a surface is discussed with respect to the distance 
from the nozzle to the impingement surface, impingement 
surface shape, shear stress on the surface and melting 
effects. Simulation of a block of ice being melted by a 
water jet is considered. The jet impinges on the ice sur- 
face, the heat transfer calculated, the change in ice sur- 
face shape due to melting is found and a new ice surface 


shape is generated. The procedure is then repeated. To 
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obtain a qualitative view of the melting process, photo- 
graphs of a block of ice were taken at various intervals 
during the melting process. Finally the combined effect 
of melting heat transfer and erosion caused by shear stress 


are discussed with respect to axisymmetric jet impingement. 
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CHAPTER I 


INTRODUCTION 


Impinging jets used as heat transfer devices, in 
general, have a great practical importance for many indus- 
trial applications. A few industrial uses are the anneal- 
ing of non-ferrous sheet metals, the tempering of glass, 
drying of textiles and paper, and turbine blade cooling. 
In fact, in any situation requiring local heating or coo?- 
ing of industrial equipment, impinging jets may have 
applications. 

Another possible use of impinging jets might be in 
the melting of ice. Mellor [42] gives some potential 
applications of ice cutting by water jets. Some of these 
are the aiding of river and lake icebreakers, deicing of 
road and runway pavements, as a cutter for the laying of 
an arctic pipeline, and as a protective cutter to prevent 
the icing of piers and pillings in harbours during winter. 
The cutting of frozen ground may be another potential 
application of impingement jet cutting. Here, both the 
mechanical action of the jet and the jet heat transfer pro- 
perties will influence the cutting process. 

There are two general cases of jet impingement. In 
the first case, a fluid jet issues into an ambient medium 


Oflassimilear fluid, liguid into liquid or gas. into gas. 
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The entrainment of the ambient fluid by the jet determines 
the flow field and consequently the heat transfer char- 
acteristics on the impingement surface. Many authors [24, 
29-38] have studied the impingement heat transfer char- 
acteristics of these jets. Some of the important effects 
are the distance from the nozzle exit to the impingement 
surface, jet Reynolds number, Prandtl] number, turbulence 
intensity and scale, and the impingement surface geometry. 
Because of the complex nature of the interactions between 
the jet, ambient fluid and impingement surface, especially 
in the region of the stagnation point, these investigations 
have been experimental in nature. 

The second case of jet impingement occurs when a 
liquid jet issues into a gaseous medium. The jet forms a 
free surface between it and the ambient medium. A diffi- 
culty in this problem is that the position of the free sur- 
face is unknown and it's location must be found as part of 
the solution. It is easier to analyze this case of jet 
flow analytically. As stagnation point fluid flow and 
heat transfer can be analyzed by boundary layer methods, 
the division of the problem into a seven al flow region 
and a boundary layer is appropriate. 

The potential flow problem can be formulated as an 
elliptic boundary value problem or as a variational state- 
ment minimizing a functional. When the potential flow is 
‘Stated as a boundary value problem the usual methods of 


solution are separation of variables and finite difference 
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numerical solutions. The methods of complex variables 

and conformal mapping may be used as a solution technique 
for two-dimensional potential flows only. As a variational 
statement, the solution may be obtained by using the 

finite element method. Solutions for the potential flow 
field may also be found experimentally. The electrical 
analogy, and the measurement of the pressure distribution 
in the flow field are two experimental methods. 

An exact solution for the potential flow of a two- 
dimensional jet issuing from a nozzle of a finite height 
and impinging on a flat surface is presented by Miyazaki 
and Silberman [23]. The solution method is a conformal 
mapping of the flow field. The velocity distribution on 
the impingement surface, important for boundary layer 
analysis, is found directly from the complex variable 
analysis. Separation of variables has been used by Scholtz 
and Trass [28] to solve the potential flow of a non- 
uniform axisymmetric jet impinging on a flat surface and 
by Sparrow and Lee [47] for a non-uniform two-dimensional 
jet. Non-uniformity exists because of the assumption of 
a fully developed jet velocity profile in the nozzle. 

The location of the free surface is found by a method of 
successive approximations. Sarpkaya and Hiriart [16] 

have used the finite element formulation of Chan and 
Larock [15] to find the potential flow field of axisymme- 
‘tric jets impinging on curved target type thrust reversers. 


Potential flow of jets impinging on curved deflectors was 
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solved by a finite difference relaxation technique by 
Schnurr, Williamson and Tatom [49]. Curved deflectors 
and thrust reversers are important in connection with STOL 
and conventional aircraft. Experimentally, Leclerc [26] 
has used the electrical analogy and Brady and Ludwig [27] 
have measured the pressure distribution in the flow field 
of an axisymmetric jet impinging on a flat surface. Once 
the pressure distribution is known, the velocity distribu- 
tion is calculated by Bernoulli's equation. 

Of the methods mentioned above, the finite element 
method seems to offer the most flexibility in solving 
these jet impingement problems. Finite elements are well 
Suited to both axisymmetric and two-dimensional flows as 
well as arbitrarily shaped boundaries. The method of 
complex variables is only useful for two-dimensional 
potential flow and not applicable to axisymmetric flow. 
Separation of variables requires simple boundary conditions 
and simple boundary geometry for calculations to be 
possible. The finite difference method requires many more 
nodes than the finite element method to adequately describe 
a curved boundary shape such as the shape of the free sur- 
face. Therefore, because of the Flexibility of the finite 
element method, the potential flow jet impingement con- 
Sidered in this thesis was solved by the finite element 
method. 

Boundary layer analysis is used to calculate the skin 


friction and heat transfer coefficients associated with a 
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particular flow. Miyazaki and Silberman [23] used the 
exact solution for the potential flow of a two-dimensional 
jet and then employed a finite difference boundary layer 
numerical technique to find the heat transfer and skin 
friction coefficients on the impingement surface. Another 
analytical method is to calculate the boundary layer para- 
meters associated with Hiemenz stagnation point flow as 
described in Schlichting [17]. This method requires know- 
ledge of the potential flow in the region of the stagnation 
point and is only valid in the immediate neighbourhood of 
the stagnation point. 

Experimentally, the heat transfer coefficients may 
be found basically by two methods. The first is by using 
a napthalene sublimation mass transfer technique. Scholtz 
and Trass [28] used this technique for the non-uniform 
axisymmetric jet impinging on a flat surface. The heat- 
mass transfer analogy is then used to determine the heat 
transfer coefficients from the measurements of mass trans- 
fer. The second method is to use a hot wire or hot film 
heat sensing device. Baines and Keffer [48] have used this 
method to determine the heat transfer and shear stress at 
a stagnation point arising from the impingement of a fully 
developed turbulent jet. 

Because of the extremely complex nature of a water 
jet impinging on an ice surface there has not been much 
analytical treatment of this subject. Yen [41] has experi- 


mentally studied a bubble-induced water jet impinging on 
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an ice surface. The effects of entrainment are important 
in this case. Yen and Zehnder [43] and Gilpin [44] have 
studied the melting of a block of ice by an impinging 
water jet. In [43] the melting effect was determined by 
the mass of ice removed by the impinging jet. In [44] two 
distinct modes of melting were noted, differentiated by 
the shape of the ice cavity the water jet produced. The 
first was a smooth cavity shape occurring at small jet 
Reynolds numbers and the second was a rough cavity shape 
occurring at larger jet Reynolds number and is characteris- 
tic of fully turbulent flow. 

Melting ice under conditions of forced convection 
has been analyzed for the case of laminar flow on a flat 
ice surface by Pozvonkov, Shurgalskii and Akselrod [19]. 
Their solution was attained by utilizing the integral 
boundary layer equations. As no pressure gradient term 
appears in the flow over a flat plate, their solution is 
not applicable to the present case of jet impingement flow. 
Other analytical treatments of forced convection melting 
of ice [50-53] also involve the case where the pressure 
gradient term has been neglected. 

In this thesis the analytical treatment of a laminar 
water jet impinging on an ice surface is discussed. 
Chapter II deals with the potential flow problem formula- 
tion, as a boundary value problem and as a variational 
statement. The equivalence between the boundary value 


problem and the variational statement is shown by the 
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methods of the calculus of variations. Chapter III pre- 
sents the finite element method as a solution to the 
potential flow problem. Chapter IV derives the boundary 
layer solution. The Karman-Pohlhausen integral method is 
used with a fourth order polynomial to approximate the 
velocity and temperature distributions in the boundary 
layer. Both the momentum integral equation and energy 
integral equation are formulated and it is seen that they 
are coupled through the boundary condition at the melting 
interface. Chapter V presents calculations based on two- 
dimensional jet impingement in the absence of melting. 
Comparisons with exact solutions and other numerical cal- 
culations is given to show confidence in the present method. 
Chapter VI deals with axisymmetric jet impingement. The 
effects of nozzle-plate spacing and impingement surface 
Shape on the surface velocity distribution are discussed. 
In addition, the shear stress on the surface resulting 
from the impinging jet is noted. Chapter VII deals with 
the melting and erosion of a frozen material. If the 
material is ice then heat transfer between the water jet 
and the ice surface is the primary ice removal mechanism. 
When the material is frozen gravel or sand, for example, 
both removal mechanisms, heat transfer and erosion will be 
important. The effect of melting on heat transfer is dis- 
cussed and an attempt is made to model the melting of a 
block of ice. Finally the combined removal mechanisms, 
heat transfer and erosion are discussed with respect to an 


axisymmetric jet impinging on a flat surface. 
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CHAPTER II 


POTENTIAL FLOW ANALYSIS 


The classical approach for the solution of many 
forced convection momentum and heat transfer problems 
begins with the computation of the inviscid flow, ignoring 
any boundary layers which may be present. The boundary 
layer is then computed assuming the flow conditions at 
the outer edge of the boundary layer are equal to the 
conditions predicted at the surface of the body by the 
inviscid flow analysis. The parameters of interest, the 
skin friction and heat transfer coefficients can then be 
found from boundary layer considerations. 

In this chapter, the steady, irrotational flow of an 
inviscid incompressible fluid is considered. The flow 
field problem is formulated as a boundary value problem 
and as a variational statement both involving the scalar 
velocity potential and the stream function. Both two- 
dimensional, using plane cartesian coordinates and axi- 
symmetrical flow using cylindrical polar coordinates 
independent of the angular coordinate are considered. The 
cartesian coordinates used are x,y and the cylindrical 
polar coordinates used are r,z. | 

The science and mathematics of inviscid flow analysis 


is very well established and the details may be found in 
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many standard textbooks, for example, Milne-Thompson [1], 
Batchelor [2], Lamb [3] and Robertson [4]. The boundary 
value problem formulation is given here for completeness 

and future reference. The variational formulation follows 
that of Kantorovich and Krylov [5] and is included to show 
the equivalence between the boundary value problem and the 
variational statement, especially with respect to the bound- 


ary conditions. 


2.1 Governing Equations 

Consider the steady, irrotational flow of an inviscid 
incompressible fluid. Let prime (') denote a dimensional 
variable and unprimed variables as non-dimensional. The 


equation of continuity “ts 


where V' is the velocity vector. This reduces to 
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for two-dimensional flow and to 
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for cylindrical axisymmetrical flow. The condition of 


irrotationality is mathematically expressed by 


Ver Xe =O Zend 
which becomes 
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for two-dimensional and axisymmetrical flow respectively. 
lo saciscy the scontinuity: equations 2.1.2 and 2.1.3, a 
stream function w~' is introduced and is defined by 
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From the condition of irrotationality, 2.1.4, the 
curl of the velocity vector is zero and from vector analysis 
the curl of the gradient of a scalar function is also zero. 
Therefore, for irrotational flow the velocity vector is the 
gradient of a scalar function. Define a scalar velocity 


potential function ¢$' as 


Vite o.v oO: ran | a 


The velocity components in terms of the potential function 


are 
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for two-dimensional and axisymmetrical flows respectively. 
Letting V be a reference velocity and D a reference length 


a set of non-dimensional variables may be defined. 
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Using non-dimensionalizations 2.1.12 the continuity equa- 
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for axisymmetrical flow. Using 2.1.12 the conditions of 


irrotationality, equations 2.1.4, 2.1.5 and 2.1.6 become 
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for two-dimensional flow and 
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for axisymmetrical flow. The velocity components in terms 
of the non-dimensional stream function yw and the non- 
dimensional potential function $ are found by substituting 


2.1.12 into equations 2.1.7 through 2.1.11. Therefore 
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in general for irrotational flows. Substitution of equa- 


tion 2.1.21 into the continuity equation 2.1.13 yields 


Vo = 0 es leed 
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for axisymmetrical flows. Substitution of the velocity 
components in terms of the stream function, 2.1.19 and 


z-!.20 into equations 2.1.17 and 247.18 fesults in 
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ine*three equations; 2.1.23, 2.1.24 and 2.1.25 are 
Laplace's equation but the equation for the axisymmetric 
stream function 2.1.26 is not. Consideration of equation 


2.4 19S yields 


Which are recognized as the Cauchy-Rieman conditions. The 
functions » and ~ are conjugate harmonic functions and two- 
dimensional problems can be solved by the methods of com- 
plex variables. 

For the axisymmetric case, consideration of equations 


2.1.20 result in 


The Cauchy-Rieman conditions are not satisfied and @# and ) 
are not conjugate harmonic functions. Although # satisfies 
Laplace's equation 2.1624, v does not, ‘equation 2.1.26. 

The problem may be formulated either in terms of the 
velocity potential by equations 2.1.23 or 2.1.24 or in terms 
of the stream function by equations 2.1.25 or 2.1.26 
depending if the flow is two-dimensional or axisymmetric. 
The velocity components are given by equations 2.1.19 and 
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2.2 Boundary Conditions 

For the derived governing equations 2.1.23 through 
2.1.26 either Dirichlet or Neumann type boundary condi- 
tions apply. At a solid boundary the velocity component 
normal to the boundary must be the same as the velocity 
of the boundary in the normal direction. Therefore if n 


is the normal unit vector 


In the two-dimensional case this results in 


Velie have tit, = Vpyn + Ve on CHeee 


When a solid boundary is stationary, Vp = 0 and using equa- 


tion-2.1.19 
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and 
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Substituting equations 2.2.4 into 2.2.3 results in 
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Equations 2.2.5 state that the tangential derivative of the 
stream function along a stationary solid wall is zero, while 
the normal derivative of the potential function is zero. 
That is, the velocity normal to the solid wall is zero. 

The boundary condition of a constant specified velo- 


city yom normal to a boundary segment is given by either 
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line through which there is no normal flow. Equations 
2.2.5 and 2.2.6 are Neumann type boundary conditions. 

Dirichlet conditions imply the stream function wv or 
the potential function » being constant along a boundary 
segment. If the boundary is a streamline the stream func- 
tion is constant along the boundary. If the potential 
function is constant along a boundary then the streamlines 
at the boundary are at right angles everywhere on the 
boundary. 

The boundary value problem for the steady irrotational 
flow of an inviscid incompressible fluid may be formulated 
by a number of boundary value problems. 

Given a flow region R bounded by a curve IT the bound- 


ary value problem for two-dimensional flow is either 
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subject to 226 


where [ = ry + Ty. For axisymmetric flow the boundary value 


problem for flow region R and bounding curved T is either 
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and again Tf = ry + ry. 

The above formulations can usually be solved directly 
by an analytical or numerical method if the bounding curve 
r is known a priori. In free surface flow problems the 
position of the free surface is not known before hand and 
must be found as part of the solution. A free surface is 
a streamline on which the pressure is constant, for example, 
a jet of water issuing from a nozzle into air forms a 
free surface. 

As the free surface is a streamline the velocity 
normal to the free surface is zero. In addition the condi- 
tion of constant pressure along the free streamline must 


nold. (oT hats is 
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on the free streamline. 
Integration of the equations of motion along a stream- 


line results in the well known Bernoulli equation 


I) 
a + + gh' = B™ eerale 


where B* is the Bernoulli constant for the streamline and 


v'2 is the square of the speed of fluid particles on the 


streamline; V'2 = Vor ty If the gravity forces can be 


neglected, assumed valid for jet flows to be considered, 
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then the condition of constant pressure along the stream- 
line implies a constant speed by Bernoulli's equation 
2.2.12. So at the free surface two conditions must be 


satisfied 


* ° 
P = pon the streamline 


and if gravity forces can be neglected the second condition 


reduces to 


V = v* on the streamline PrP T4 


2.3 Variational Formulation 

There are many methods available for the direct solu- 
tion of the boundary value problems given by equations 
coce/* through =2,2:10."*O0f*these are the exact solutions of 
complex variables for two-dimensional flows and separation 
of variables methods for either Poe dinensional or axi- 
symmetric flows. Some approximate methods are the finite 
difference relaxation techniques and the electrical analogy 
graphical method. 

Another method of solution would be to try and find 
‘an equivalent form of the boundary value problem. As is 


well known from structural and continuum mechanics, varia- 
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tional methods using the calculus of variations offer an 
alternative formulation for boundary value problems [6]. 
One method of solution of the variational method is the 
Finite Element Method [7,8,9]. Zienkiewicz and Cheung 
[10] first applied the Finite Element Method to the solution 
of field problems as an extension of it's use in solving 
problems of structural and continuum mechanics. Since then 
the Finite Element Method has been used to solve many 
different problems of fluid mechanics including the flow 
of an ideal fluid, for example [11,12,13,14]. Reference 
[8] contains an extensive recent list of finite element 
solutions to fluid mechanics problems. It will be instruc- 
tive to show the equivalence between one of the boundary 
value problems stated and its associated variational state- 
ment especially with respect to the boundary conditions. 
The derivation follows that of Kantorovich and Krylov [5]. 
Consider the two-dimensional boundary value problem 
formulated in terms of the potential function, equation 


2.2.7. The variational formulation considers the integral 
ee ey eer + Ear, dxdy - p f ov"ds ee 
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The boundary conditions are exactly the same as those used 
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in the boundary value problem, i.e., ¢ = 6” on ry and se = 
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boundary conditions which make the integral I a minimum. 
menetore | is a function of @ or I = I(s). Let o(x,y) 

be the function which minimizes the Bntegral 1. Now, let 
the set of all other admissible functions satisfying the 


boundary conditions, $(x,y), be defined by 


o(xsy) = o(x,y) + en(x,y) ee 


where € is a parameter and n(x,y) is arbitrary. Consider 
the Dirichlet boundary condition, » = * on Ty. Since both 
@ and ¢ must satisfy this boundary condition equatvon 2.3.2 


results in 
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then obtained by evaluating 


al(¢ + en) c 
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and using Green's Theorem to change an area integral to a 


line integral, equation 2.3.6 becomes 
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n(x,y) = 0 using equation 2.3.3. Therefore, the line 
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may be combined. So 


2 2 
a1 , 
ae = - ff n(2$+ 28yaxay 
€=0 R OX oy 
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As the area integral and line integral are independent and 


that n(x,y) is arbitrary, the result is 
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and the boundary value problem is recovered. A similar 
Procedure can be used for the other cases. The integrals 


are listed below. 


Two Dimensional-Stream Function 


I(y) = 5 i (34)% 4 oe dxdy - 9 I vids 7.2.3.9 
| 2 


where 


Axisymmetric Potential Function 


1(6) = pm it (22)% 4 (2) “\ rdrdz = 2pm i ov"ds 2.3.10 


with the same boundary conditions as in the two-dimensional 


case. 


Axisymmetric Stream Function 
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with boundary conditions 


* 
y= wv on ry 
lesvrings 
ron vVoon ly 


EGUatTONS? 2.3.4 Ae 2.3.9 are equally useful for two- 
dimensional flow and the choice of which one is most suited 
to a particular problem depends on the type of boundary 
conditions that are present. However, equation 2.3.11 
wsenOteas Useful as 253510 due’ to the presence of the 
radial coordinate r in the denominator of the first integral 
of the equation. To evaluate this integral it would be 
necessary to use numerical integration which is time consum- 
ing and susceptible to errors. Equation 2.3.10 will be 


used for all axisymmetric flow problems. 
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CHAPTER III 


THE FINITE ELEMENT METHOD 


3.1 Background 

The development of finite element techniques origin- 
ated from classical approaches to structural analysis. 

The sciences of structural and continuum mechanics have 
used finite element techniques to solve many complex pro- 
blems. Recently, the finite element method has been used 
to solve many field problems including those of potential 
flow in fluid mechanics [7-14]. 

The finite element method has been used extensively 
in recent years because it has, in general, several out- 
standing features. Listed by Chan [11] they are: 

- Non-homogeneous and anisotropic problems are 
easily treated. 

- Elements can be graded in size and shape to follow 
boundaries of arbitrary configuration. 

- Once a computer program has been developed, pro- 
blems of a similar nature can easily be solved by supply- 
ing the appropriate new input data. 

The general method of solution of fluid flow problems 
using the finite element method is as follows: 

- The entire flow region is divided into a set of 


elements, interconnected at various nodal points. A number- 
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ing system is chosen to assign values to the nodes and 
elements. 

- An element interpolation function is chosen specify- 
ing the variation of the field variable within each element. 
In the cases considered the field variable is either the 
Stream function or the potential function. 

- The contributions of each element to the total flow 
pattern are assembled by a Ritz technique resubting in’ a 
System of symmetric, banded, linear simultaneous equations. 
The number of equations in the system equals the number of 
nodal points in the element mesh and the solution of the 
System gives the value of the field variable at each nodal 
point. 

- All related physical properties are evaluated from 
the nodal values. 

The Ritz technique is a procedure for changing a 
continuous medium problem into an approximate lumped para- 


meter system. Consider the integral 
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over the domain R. The object is to minimize I(u); that 
is find the function u from a set of admissible functions 
which makes the integral a minimum. This is accomplished 


by selecting an appropriate trial family of solutions Uns 


where U, is given by 
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The u. are the undetermined parameters and the Y; are the 
coordinate functions or trial functions. The coordinate 
functions are known before hand but the u. are to be deter- 
mined so to make UL, the function which minimizes I(u). 


The minimization is achieved by finding the quantity 


=70 Seale 


The result is a set of linear equations which enables the 
best approximation to the true solution from the family of 


trial solutions to be found. 


3.2 The Finite Element Method Applied to Jet Impingement 


(a) Finite Element Formulation 

The following is a description of the finite element 
method as applied to a jet impingement problem. The formu- 
lation follows that of Chan [11] and Chan and Larock [15]. 
Sarpkaya and Hiriart [16] have used this formulation to 
salve the jet flow associated with aerodynamic thrust 
reversers, a very similar problem. 

Consider the flow of an axisymmetric jet impinging 
on a flat plate. The flow field is shown schematically in 


Figure 3.1. As the flow is axisymmetrical the potential 
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function formulation must be used and the appropriate 
integral is equation 2.3.10. The procedure followed here 
is identical with velocity potential and stream function 
formulations in two-dimensional flow. Since the flow is 
Symmetrical about the stagnation streamline AB, only one- 
half of the flow region need be considered. 

First the flow region is divided into a triangular 
element mesh, and a numbering system assigned to the nodes 
and elements. The pattern chosen is shown in PaLgure 3.2. 
There are a number of reasons for a pattern of this type. 
A sequence is established for the numbering of both the 
elements and the nodes. The elements shown in the figure 
are in blocks of six. All of the blocks of six elements 
are identical in the flow field. This enables the flow 
field to be modeled by as many elements as desired PUSt DY 
adding more blocks of six elements. The node numbering 
System also follows a distinct pattern. It is possible to 
generate the numbering system for the whole nodal field by 
Starting with the node numbers of the first six elements. 
Therefore, a computer program may be written to generate 
the “TR system for the flow field just by inputting 
the numbering pattern for the first six elements and the 
number of blocks of six elements desired. 

A six node triangular element is used as the basic 
element of the flow field. This enables a quadratic varia- 
tion of the field variable, the potential function $6, within 


the element. Introduction of triangular area or natural 
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coordinates leads to simplification of the algebra 


involved. The area coordinates ar are defined by 


where a is the area of the triangular element and the A. 
are. the areas of the subtriangles as shown in Figure 3.3. 
Each element has three corner nodes (numbers 1, 2 and 3) 

and three mid-side nodes (numbers 4, 5 and 6). For example, 
the side connecting node 1 to node 2 is described by 53 = 0 
with Ey varying from 1 to O and Eo Varying from QO to 1. 
Obviously, a constraint on the area coordinates defined by 


equation 3.2.1 75 
oy + So hi 63 = ] Oe cenc 
The quadratic variation of the potential, go”, in ele- 


ment m is given in terms of the potentials o; at the node 


points. 
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The derivation of the element interpolation function X, 

is given by Huebner [8]. The original polar coorindates 
htezdnof a point in the triangle are linearly related to 
the area coordinates (E,; Eos E3) by the system of equa- 


tions represented in matrix form by 


zy Z» Z. Ey Z 
ry r5 r3 bE = Yr Sige nO 
] ] 1 b3 1 


where (r,, Zz.) are the coordinates of node i and i=(1,2,3). 


Solving equations 3.2.5 for (E); Ep: E3) results in 
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j = (2,3,1); k = (3,1,2). In indical notation equation 3.2.6 
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As the variation of gm is quadratic in the element, 
the variation of the velocity will be linear within the 
element. The velocity components are given by equation 
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The boundary conditions will be considered later. To mini- 
mize 3.2.10, differentiation with respect to the nodal 


Values oF as discussed in equation 3.1.3 is required. 
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This expression gives the contribution of the mth element 
to the overall flow field calculation. The integrals are 
evaluated in terms of area coordinates and transformations 


given by equation 3.2.5 are used. In short hand form 
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andwies Lhtolésagss TAtos6: 

As the integration is carried out in area coordinates, 
which are the same for each element, the integrations need 
only be done once. The element matrices SAG; and SLAS LON 
axisymmetrical elements and Si and Site for two-dimensional 
elements are listed in Appendix I. Also given is a sample 
integration. 

Now, essentially a bookkeeping operation is required 
to add up the contributions of each element to form a global 
set of simultaneous equations. For the minimum by equation 


3.1.3 
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Note that the factor 2pm is common to all the terms and may 


be divided out. The resulting matrix equation is of the form 
[SAI{¢} - [SLA] = 0 Sie". als 


In terms of structural analysis the matrix [SA] is the stiff- 
ness matrix and [SLA] the load matrix. This set Of equa- 
tone is symmetric,banded and linear and can be solved 
efficiently by a direct Gaussian elimination technique. 

A saving is realized if the bandwidth is minimized. 
‘The half-bandwidth is determined by the node numbers of a 


given element. It is equal to one, plus the difference 
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between the largest and smallest node numbers in the ele- 
ment. As shown in Figure 3.2 the numbering of the nodes 
is such that the minimum bandwidth is realized. The half- 
bandwidth for those elements is seventeen. Because the 
matrix [SA] is banded and symmetric its size is [NN x NBW] 
Where NN is the number of nodes and NBW is seventeen, the 


half-bandwidth. 


(b) Boundary Conditions 

The boundary conditions appear in the matrix [SLA]. 
Referring to Figure 3.1 they are as follows. 

(1) Symmetry Line AB: This is the stagnation stream- 


line and there is no flow normal to this line. 
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(2) Impingement Surface BC: There is no flow normal 


to the solid wall. 
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(3) Outlet Section CD: Assume that there is uniform 
flow at the outlet section. The velocity in the z direction 
is zero and the velocity in the r direction is assumed 


constant. 
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Veseeere She toRo nas YF is 3.2.18 


(4) Free Streamline DE: As this is a streamline 


there is no flow normal to it. 
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The other condition of constant velocity on the free stream- 
line is arrived at through an iterative procedure. An 
initial guess as to the position of the free streamline is 
made and the potential flow equations are solved. The 
velocities on the free streamline are then calculated and 
the position of the free streamline is adjusted to give a 
better approximation of the free surface which satisfies 
the constant velocity boundary condition. A more detailed 
discussion on the adjustment of the free Surface is given 
in the next section. 

(5) Nozzle Streamline EF: Again there is no flow 


normal to a solid wall’ and 
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(6) Inlet Section FA: Assume uniform flow at the 
‘inlet section. The velocity in the r direction is zero and 


the velocity in the z direction is constant. 
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The problem as formulated for jet impingement on a flat 


plate is 


Sacaize 


Equations 3.2.22 constitute a Neumann boundary value problem 
Which has a unique solution to within an arbitrary constant. 
To give a totally unique solution a value of the velocity 
potential may be specified at a point. Consider the bound- 


ary condition at the inlet section FA. 
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Integrate 3.2.23 to get 
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Where F(r) is some function of r. From 3.2.24, F(r) can 


only be a constant at most on the inlet section. At the 
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inlet section z is also a constant. Therefore an alterna- 


tive boundary condition for uniform flow at the inlet is 


do = C2 on FA BA1Z, 25 


As this is a Dirichlet boundary condition the inlet section 
of the boundary IT becomes r,. The only non-zero part of 
ry is the outlet section CD. The constant C3 is arbitrary 


and may be set equal to zero. 


(c) Finding the Free Surface 

The position of the free surface is not known a priori 
and must be found as part of the solution. If gravity 
effects are neglected the condition of constant velocity 
is satisfied when the position of the free surface is found. 
An iterative technique, following that in [16] is used. 
An initial guess is made for the position of the free sur- 
face. The potential flow finite element equations are 
solved and the velocity at the nodes on the free surface 
determined. The velocity is compared with the velocity 
calculated for the node at the lip of the nozzle Ves The 
lip of the nozzle is point E in Figure 3.1. The velocity 
on the streamline should be constant and equal to Ves Nay 
the velocity calculated at a node on the free surface is 
larger than Ves the node is moved in a direction normal to 
the streamline so as to decrease the velocity on the next 


iteration. If the calculated velocity is smaller than Ve 
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then the node is moved in a direction so as to increase the 
velocity on the next iteration. Four cases of node move- 
ment are shown in Figure 3.4. In each case the node N is 
moved to the position N'. The movement is based on general 
continuity requirements. The distance moved normal to the 


Streamlines is given by 


V 
An = + mY oPh B26 


where. is a parameter similar to a relaxation parameter. 
The square of the ratio of the calculated nodal velocity 
Vine to the velocity at the lip of the nozzle Ve is used to 
accentuate the velocity difference. 

This velocity comparison and Subsequent movement of 
the free surface node is done for all nodes on the free sur- 
face. Once a new position of free surface is found, the 
finite element solution procedure is repeated, the velocity 
on the free surface again found and the velocities compared 


to v The iteration is stopped when the calculated velo- 


FE: 
city agrees with the actual velocity to within a certain 
percentage. Usually the agreement limit was set between 
One and two percent. 

As pointed out by Sarpkaya and Hiriart [16] converg- 
ence to the proper solution can be ensured with an appropri- 


‘ate choice of the parameter A. They suggest a value of 


0.015. This value was used and convergence always occurred, 
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The most important information derived from the 
potential flow analysis is the distribution of velocity 
along the solid wall, segment BC in Figure 3.1. This 
distribution is directly related to the pressure gradient 
in the boundary layer which grows from the Stagnation point. 
The velocity distribution is required before the solution 


of the boundary layer may be considered. 
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CHAPTER IV 


BOUNDARY LAYER SOLUTION 


4.1 Governing Equations 

Consider the laminar boundary layer growing from a 
stagnation point of an incompressible warm water jet 
impinging on an ice surface. A schematic of the flow is 
shown in Figure 4.1. Because the POS@tion of the ice- 
water interface is not fixed but dependent on time and the 
rate of heat transfer, the problem is an unsteady one. 
This complication is removed by considering a coordinate 
System which moves with the melting ice front. The problem 
is now assumed to be a quasi-steady state problem. 

Neglecting viscous dissipation and assuming that the 


radius of curvature R' of the ice surface is sufficiently 


6 
large such that the condition an Cris Satisfied, [18] 
the boundary layer equations for two-dimensional flow 
(n = 0) and axisymmetrical flow (n = 1) of a constant pro- 


perty fluid are 
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Figure 4.1 Boundary layer on an ice surface 
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MOMENTUM: 
t i] i] 2 1 
- = feu ee = iy" qv feats 0 5 A. 1. 
aX ay dx dy 
ENERGY: 
iy) La ee ea re 
9x" ay. ay 2 ol. 
The boundary conditions are 
at y= 0 flow conditions u'n=10; v's F 
thermal conditions T = To qs 
' Zot 
BENT yHUs 5 flow conditions ous =0; 4 ie =(0; ou = 
dy oy 
i] 
at y' = 6! thermal conditions 2 = 0; ie SRG. Wad 
t ay! dy 2 Toe) 


Another boundary condition is realized at the melting sur- 
face. Assuming that the ice is at a constant temperature, 
the melting temperature 0°C, an energy balance at the melt- 


ing interface results in 


pelv, =k 2u] , 4.1.5 
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where Vo is the injection velocity due to the melting ice. 


The following non-dimensional variables are introduced: 


5, T-T 
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Where V is the velocity of the impinging jet and D is either 
the jet diameter or slot width depending if the flow is axi- 
Symmetrical or two-dimensional. With these substitutions 


the boundary layer equations become 
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and the boundary conditions become 


at’ ys =20¢tue =f Oe byt= 4 


0 
0 = 0 
ou at 
aty=6 3 -——= 0; =.= u = U(x) 4.1.10 
u oy 2 
oy 
at y 26 2S 20% Bao! AE 5 o= 1 
t ay ay 
The energy balance at the interface, 4.1.5, becomes 
Be as eee) 
Mo # Provoy y=0 aol 
Where Ste is the Stefan number defined by 
Ta To 
Ste = SM xe renee eal Reg Ie 


The Stefan number gives the ratio of specific heat of the 
water to the latent heat for the melting phenomena. It is 
always positive (tT. > Ty) for melting. As the temperature 
range for the water jet is 0°C to 100°C the Stefan number 
is limited to the range 0 < Ste < 1.25. 

| The method of solution of equations 4.1.7, 4.1.8 and 
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4.1.9 with boundary conditions 4.1.10 and Adee Peis ethe 
Karman-Pohlhausen integral method [17,18]. A fourth order 
polynomial is assumed for both the velocity and temperature 


profiles in the boundary layer. 


4.2 Integration of the Boundary Layer Equations 
The boundary layer equations are going to be inte- 


grated with respect to y from Ya= 0) tO rye = 26, 
CONTINUITY 


e  (r"u)dy + ie =< (r"v)dy = 0 
9 ox 9 oY 


Performing the integration and applying the limits 
] 0) n 
= a F Ge (r-u)dy a2. 


Both 6 = 5 and 6 = 6, will be considered when using 


the integrated fOFVMVon the continuity equation 402 21- 
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2 2 n 
ou CU oOMU ee UN cy eu ser kyl ny? d(uv) 
oY y 7 ae ay Hr GX) hh ox ar ae aes ay 
ue dr" 
using the continuity equation 4.1.7. The term Batdsce 
r 
vanishes for the two-dimensional case (n = 0).. The first 
integral then becomes 
pu +S (mu? yay + pu otuv) gy = - ued ro) 
jr OX 0 oy lt dx 2 
6 2 
dU u u u 
aCe veh (ota ay 2 UNE 
0 U 
Where boundary conditions 4.1.10 were used and 55 is the 
momentum thickness defined by 
Sy u 
65 a f U (1 = y) dy Atk 
The third integral of equation 4.2.2 evaluates to 
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using boundary conditions 4.1.10. Combine the second 
integral of 4.2.2 with the left hand side and multiply by 


-1 and the momentum equation becomes 
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The final form of the integrated momentum equation is 
arrived at by dividing by U and using equation 4.1.11 for 


Vo to obtain 
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Using the continuity equation 4.1.7 and boundary conditions 


4.1.10 the left hand side integrates to 
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where 64 is an energy thickness defined by 
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The right hand side becomes 
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Multiplying by -1 and substituting for Vo from equation 


4.1.11, the integrated energy equation is 
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The effect of melting is seen to appear in both equations 
4.2.5 and 4.2.7 with the presence of the Stefan number. If 
the Stefan number is zero these equations reduce to those 
given in [17,18] as the integrated boundary layer equations 
of an incompressible fluid in the presence of a pressure 


gradient for either axisymmetric or two-dimensional flow. 


4.3 Velocity and Temperature Distributions 
A fourth order polynomial is chosen to represent the 


velocity and temperature distributions in the boundary layer. 


(a) Velocity Profile 


Assume the form of the velocity profile as 
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The constants ay to a, are evaluated using boundary condi- 


tions 4.1.10 and evaluating the momentum equation 4.1.8 at 


y= 0. 


te 


TSA 


2nortaups ntod of reagee os nae2 ae ve doa 
TT .aeden netat? oas to. slain Real Sat 888 
seodd ot g0tbax snotiquee aepit over 2f otimise 
enorseves +syel uisboued basexeadat afd 28 [81007 pa 
otezonc & to soaeest@ ome ai bret ers a 
wort reno Kenan 0 ais singe tae sadsis 40% brbe 
“ODE OTRS > 

enotsudrnae ht meron on bos xateotov 

sit dng2zs71q9" oF aszors et Tetmone tog +abvo ddayot A 
yet yrohodod edt a apr reer = use 
ce x 

sorter’ ystooty (oy 

26 otto yo hooTew ads, 0) sah ons, smuaaA 

aS a 


a) ae (Tig Teli + Tiss * pore - ou =u 


a 


“bie wrsbrued gnigu betsuleve ove y5 02 yo esentenoa aAT 
$h.BH.b AotPeMPS munamoM att PotIeuleve brs (OTS fh id 
: ; ; ¥ 7 . ea | 


where 


Ss 


at 


at 


at 


a) ay + 2a + 3a 3 + 4a 4 


2 
dU 
¥en Le = 
u ] aye 


2 . au 
yu2 OF nN 0; Yo y 
sO 6A ay - 2a, = oh 

' Yo%u 


0 


98 


<3 al) 


Cb) 


43.6) 


aod) 


«3(e) 


os 4 


(aje.e.6 


(dye 6.8 


(OE cE 4 


eos eh eye eae ‘ 


4 Le at ele sacitie 


er ere ae 


99 


and 


ay aS. 


The five conditions 4.3.3(a) to (e) are sufficient to 
evaluate the five unknown constants ay to ay. The system 


of equations may be represented in matrix form as 


] 0 0 0 0 ay 0 
0 ] 2 3 4 ay 0 
0 0 ] 3 6 ao = 0 4:53 
0 1 T ] ] a4 ] 
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where 


= 3 4 
F(n,) a ény ri 2nj RE ny 
a 2 3 4 
G(n,) = ny - 3nq + 3nj_- n, 4.3.8 


3 
H(n,) = én‘ - 8n) + any 


The functions F, G, and H are shown in Figure 4.2. 
Velocity profiles for a range of A, the pressure gradient 
parameter, with of = 0 are shown in Figure 4.3 and the 
profiles for a range of ho when A = 6 are shown in Figure 
4.4. The effect of melting on the velocity profile is very 
small as indicated by Figure 4.4, The parameter X1, is a 
melting parameter and is related to the Stefan number. 

For the case of no melting, ho = 0 these profiles are 
only valid in the range, -12 si plc il 7 1S. iiemlower 
limit is the point of separation in this theory and the 
upper limit corresponds to the case when the velocity pro- 
Files pop [16]; that. 1s T > 1 for a value of ny between 0 
and 1. The condition, 7 > 1 is physically unrealistic for 
steady isothermal flow. When melting is present this condi- 
tion is modified slightly to -12 < A < 12(1 + 29). Also, 
for the case of no melting the velocity profile reduces to 
the one presented in [17] in his discussion of the Karman- 


Pohlhausen integral method. The case of melting on a flat 
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Figure 4.2 Velocity and temperature profile functions 
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Figure 4.3 Velocity profiles , Ste. = 0 , effect 
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Figure 4.4 Velocity profiles , A= 6 » effect of melting 
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plate was studied by Pozvonkov, Shurgalskii and Akselrod 
[19]. The above profile reduces to their profile when 


A = 0, that is for the case of zero pressure gradient. 


(b) Temperature Profile 
A fourth order polynomial is also chosen to represent 


the temperature profile. Assume 
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To evaluate the constants, boundary conditions 4.1.10 


and the energy equation 4.1.9 evaluated at y = 0 are used 
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2 
at y= 643 pe ae a5 = 0 
oy 
SO bs = 2 3b. + 6by = 0 AISLUT Cc ) 


at ye 643 n= ts © =] 


SO by + bo + b3 + Dy = |] 435d) 
2 
00 ih. aoe 
at y= 03 9 W!=005 bv Hg tetra Pes 
0 dy|y=0 We ay? y=0 
SO 6A_Prb, - 2b. = 0 A311 (2) 
where 
V,6 
i Oa 
dg = 6 4 312 


This system of equations may be represented in matrix 


form as 


] 0 0 0 0 bo 0 
0 ] 2 3 4 by 0 
0 0 ] 3 6 b, = 0 i aS 
0 ] ] ] ] b. ] 
0 6AQPr -2 0 0 by 0 
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Comparison of the system 4.3.13 with the system 4.3.6 
reveals that replacing ny by ns do by AgPr and setting A = 0 
makes the two systems identical. Therefore the temperature 


profile is given by 


0 = =a -Pr {F(n) + AgPrH(n) } 4.3.14 


where the functions F and H are those given by equation 
4.3.8. Figure 4.5 shows the temperature profiles for 
various Stefan numbers. The parameters ho and MG are melt- 
ing parameters related to the Stefan number. From equation 
4.1.11 for the melting velocity, evaluating the derivative 


from equation 4.3.14 and using equation 4.3.12 one obtains 
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and solving for AgPr 
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Figure 4.5 Temperature profiles , effect of melting. 
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4.3.4 and 4.3.40 


9 ae: Bias 1 
ho = 7 o = MA; rn Wy 
where 
5 
m = = 4.3.18 
u 


is the ratio of boundary layer thicknesses. 
-With the velocity and temperature profiles 4.3.7 and 
4.3.14 the thickness parameters defined by equations 4.2.4, 


4.2.3 and 4.2.6 may be evaluated. As oF, and 6, are not known 


E 
51 6 53 
only the ratio's <r oe and Z may be found. Define 
Vv) u t 
6 ] 
] u 
A, =x =f (1 - a)dn 
] oF 0 U ] 
6 ] 
ieee uyu 
u 0 ' 
6 1 
3 u 
A, =x =f 7 (1 - @)dn 
3 Oy 0 U 
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i0 1ieo US qos .c0 


ee — 1, {3 Boe “sre Mo 
361.4 
156 .! 0 
+ 722 Ag + | 4.3.21 


For the case of no melting, do = 0, these expressions reduce 
to those reported in Schlichting [17]. To evaluate A, a 
relationship between n and ny must be established. The 


relationship is 
ny = mm Aid wee 


from equations 4.3.2 and 4.3.10 where m is the boundary 


layer thickness ratio in equation 4.3.18. The integral A 


can now be evaluated as 


] 


phe, Boundary Wver this 
3 (1 + so eee) 


{C,(m) + AC, (m) + AgC4(m) 


A_PrC,(m) * AgPrAc,(m), 4 ApPragc,(m)} eee 
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C, (m) 7 


| 
S| 


m 
Co(m) = 99 - a4 * Tes0 ™ - a6” 


Ce Gu = ne : an” : mn 4.3.24 
Cy(m) = 5 - 7 : 5 
Cg(m) = 35m - ip ‘ 


This value of A. reduces to the one reported by 
Schlichting [17] for a case of no melting, Ag = Ag = 0, 
on a flat plate, A = 0, for a Prandtl number greater than 
one. When the Prandtl number is greater than one the ratio 


of the boundary layer thicknesses is less than one. 


4.4 Solution of the Integral Equations 

The solution procedure described below is that out- 
lined in Schlichting [17] extended to solve both the momen- 
tum and energy equations. First, the solution procedure 


for the momentum equation will be described. 


KS.E.8 


oO ‘ + gh i $ a 
neds -~vedwiong ssanun ack . no ie 
orien odd sno medt teteeng et edwin Fabasxd ‘a on 


oo sis aot * seatqniotsd ape woes is 


¢ 
i ai 


ie | ano tau pl pire ws eity! Fo notiutog vi 
ine ana: at witsd beatroeed “wusbdsong aiguten® att | 
| ~nsinan ont? God evfos oF bsbastxe [xt] gntsaor gio’, at bonth 

-subssena notsutoe 9d .Sent9 -2aoktoupe vpreRe bas. mus 
sil al ad thw idee thi marae aes 7 


Li A a 
: ore! i 2 j ry ' 7 @ - “se 
~~" ae ; re SP td >. : 


7] 


From the velocity profile 4.3.7 and the temperature 
profile 4.3.14, the velocity gradient and temperature 


gradient at the wall can be evaluated 


A 
ou ae U2 aoe, 4.4 
oy|y=0 Sy p ge ae 
0 
00 “ 2 
ayly-0 ~ (1 + Pray) a 


Substituting those expressions into the momentum integral 


6 A 
equation 4.2.5 and noting that ~ a =" results in 
é 2 
n A 
LE piu Be eee (eles aly = 1) ae eo bee. 
yn dx 2 dx Ay 6 ‘eta A Pro, (1 + Prig) 
2 
Multiplying by 6. and noting A, = =< 
u 2 Oy 
2 A 
uy Dy Ue oh, CU ee nai Se Se 2 
dx Pear meni ae Pt). 
2A or" dx 2 A, dx A, 1+ a0 Brite at) PAG 
Introduce the substitutions 
ar 
Z = 65 4.4 


wae ear ia dU 
Kas a> eq tt 9 ay 4.4. 
The function k is a known function as 
52 
pM SA LU eV 
k = ea OF ax AA 4.4, 
u 
from equations 4.2.3 and 4.3.5. Multiply the momentum 
equation by Ay and expand the derivative to get 
A 
dzeg2eh, Week Ube 2ste 
dx? Tey Prm(? + Pri.) 
1 + Ag 0 
\ dn" 
Sheek ear ale tela, LEE 4. 
Ay dU yn 
dx 
The last term in { } disappears for the two-dimensional 
case. An alternative form for this last term is 
drt ; 
U dx _#Uz tdr 
du a k = vn dx 4.4, 
dx 


n 
by equation 4.4.4. To eliminate the derivative, i, a 


simple transformation is used. Let 
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z* = z(r")* 4.4.8 


and the momentum equation 4.4.6 becomes 


A 
2 Le ae A 
dz 2(r") 6 1 
e et ——— + 6) ) - k(2 + —) 4.4.9 
dx U 2 7 fey 0 Ao 


using equation 4.3.15 for the Stefan number and equation 
4.3.17. Note that for the two-dimensional case n = 0, 
Z = 2” and equations 4.4.6 and 4.4.9 are identical. 

As the integral energy equation 4.2.7 is very similar 
in form to the integral momentum equation 4.2.5 the same 
procedure is used to obtain an equation similar to equation 


4.4.6. The transformed energy equation is 


dWact2 (1 + Ste) 2 ~ ind iedr vay 
dx U pag(llpsted (+ Prag)? ~ } * au ax | ue ak 
dx rn 
where 
ae ts 
W= 63 4.4.11 
T= 92 dus y au 
ee 4.4.12 


Again to remove the derivative of (r"), the transformation 
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similar to the one used for the momentum equation is intro- 


duced. Let 


w* = Wr") 4.4.13 


and using equation 4.3.15 for the Stefan number the integral 


energy equation becomes 


dw* 2(rry4 


ns 2 1 
ave 2 8a Fea + 64) - 7 4.4.14 


A Runge-Kutta technique [20] is used to solve equa- 
tions 4.4.9 and 4.4.14. The right hand side of both equa- 
tions are functions of the pressure gradient parameter A, 
and the boundary layer thickness ratio m for a given Stefan 
number, velocity distribution and surface shape. In order 
to start the integration procedure the initial values of A 
and m must be known at the stagnation point. Equations 
4.4.6 and 4.4.10 exhibit a singular behavior at the stagna- 
tion point where the velocity U is zero. As the derivatives, 
2 and a are finite at the stagnation point, the terms in 
the brackets { } on the right hand side of each equation 


must be zero at the stagnation point. Therefore two stagna- 


tion conditions exist and are given by 
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and 


dr" 
Ay fo a XgPr + sao - L-L au an = 0 4.4.16 
dx 


For the two-dimensional case, n = 0, the terms involving rn 


n 
and ae vanish. For the axisymmetric case, n = 1, and this 
term must be evaluated at the stagnation point, x = 0. 


Assume that 


dr" 
lim #v Me 4 
io Ondu An ee ] he a 
dx 


In the region of the stagnation point the velocity may be 
represented by U = 8x where 8 is a constant. Also r ~ x 
near the stagnation point. With these approximations the 
assumption of equation 4.4.17 is true. The stagnation con- 


ditions now become 
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These two equations are solved simultaneously to 
find the starting values of the parameters A and m. The 
values determined depend on the Stefan number and whether 
the flow is axisymmetrical or two-dimensional. Figure 4.6 
shows the variation of A and m with Stefan number for two- 
dimensional flow. Figure 4.7 shows the same variation for 
the axisymmetrical case. Schlichting [17] reports the values 
of A at the stagnation point for the case of no melting and 
they coincide with the values found from the solution of 
equation 4.4.18. 

Once the values of A and m are known at the stagna- 
tion point, call them No and My> a Runge-Kutta integration 
technique may be used as follows. The first assumption is 
that the boundary layer thickness ratio Mo is constant [19]. 
This assumption was found to be acceptable as the functions 
involved in the boundary layer equations are only weak 
functions of m. The integration of the momentum equation 


4.4.9 is carried out step by step in the form 


dz* 


ee ee Pa 
254A Z5 + (|, a an 4.4.19 


ae Z 


This is a Runge-Kutta method of order two-[20]. To start 
a ad dz* 
the iteration values of Zo3 2 and I, must be known 


As Ag and My are known the value of k at the stagnation 


point, Ko» may be found from equation 4.4.5. Then Zo is 
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Figure 4.6 Two-dimensional starting values 
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Figure 4.7 Axisymmetric starting values 


78 


0.55 


0.5 


0.45 


79 


given by equation 4.4.4 and “a3 by 4.4.8. For the axi- 
symmetrical case zy is zero as r" is zero at the stagnation 


point. In the two-dimensional case 25 = Zo: To evaluate 


a 
oe 0 equation 4.4.8 is differentiated and evaluated at the 


stagnation point 


For the axisymmetrical case, n = 1, and value of 


r" at the stagnation point is zero so 


dx lo 7 9 axisymmetric case 4.4.20 


For the two-dimensional case n = O and 


The general form of the differential equation for the two- 


dimensional case is 
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not known explicitly but is known implicitly through the 
pressure gradient parameter A. To find & 0° EL top tals 


rule is used because of the singularity at the stagnation 


point.— So 


dF dy 
| f dk dx? 
ax|0° ~ dE yyy 
(le-sae) Gx ado 


For the stagnation points considered in this study 

the velocity distributions in the vicinity of the stagna- 
2 

tion point are represented by U = 8x. Therefore ot = 0 


dx 
and 


= 0 for two-dimensional flows 
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The linear velocity distribution near the stagnation point 


is valid for blunt bodies only. 
dz” n 
Using Ng and ™) dx}7 is calculated with r F and ul, 


by equation 4.4.9. Now Z3 can be calculated from equation 
4.4.19. Using equation 4.4.8, Z, can be found and from 
equation 4.4.4 Ky is found. As k is a universal function 


of A by equation 4.4.5 the equation 
K(A) = k, = 0 4.4.23 


can be solved to give the pressure gradient parameter at 
the new position 1. Once Zz, has been found all of the 
relevant parameters, the boundary layer thicknesses can 
be found using equations 4.4.3 and 4.3.19. With a new 
value for the pressure gradient parameter the procedure is 
repeated until the whole boundary layer has been calculated. 
A similar technique is used to solve the energy equa- 
tion and the result is the thermal boundary thickness Sy. 
The assumption of constant boundary layer thickness ratio 
is checked as follows. In the Runge-Kutta solutions the 
primary output variables are z and W. These were deter- 
mined independently of each other. They are related as 


follows. 
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and solving for m 
IN 
m= (5 Hl? 4.4.24 
z 
oe 


This value is then compared with the value found at the 
stagnation point to check whether the boundary thickness 
ratio was constant. 

Appendix 2 gives a description of a complete solution 
procedure, potential flow plus boundary layer calculation, 


and also the computer program used for the calculation. 
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CHAPTER V 


TWO-DIMENSIONAL JET IMPINGEMENT ON 
Rerun. SURFACE 


For the case of two-dimensional flow when a free 
surface forms, an exact solution for the inviscid incom- 
pressible jet impinging on a flat surface is available by 
using the complex variable and conformal mapping techniques 
given by Michell [21] and Ehrich [22]. Miyazaki and 
Silberman [23] used the results of [22] for the potential 
flow field analysis of the impinging two-dimensional jet 
and then used a finite difference technique to solve for 
the skin friction and heat transfer in the laminar boundary 
layer region near the stagnation point. The results of 
[23] will be used to verify the validity of the finite ele- 
ment potential flow solution of Chapter III and the integral 
boundary layer solution of Chapter IV. A discussion of the 
methods of solution of [23] and comparison with the finite 


element method and integral boundary layer solution follows. 


5.1. Exact Solution of the Potential Flow Field 

The methods of complex variables and conformal mapp- 
ings may be used to determine the exact flow field solution 
of an incompressible, inviscid two-dimensional jet imping- 


ing on a flat surface [21,22]. A schematic of the flow 
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field using the non-dimensionalization of equations 2.1.12 
is shown in Figure 5.1. The relevant results for use in 
a boundary layer analysis are reported in [23] and are 


repeated here for future reference. They are 
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where H is the dimensionless nozzle-plate spacing, U. is 
the velocity on the impingement surface or mainstream velo- 
city for the boundary layer, Ue is the velocity on the free 
surface and y is a parameter defined by equation 5.1.2. 

For a given nozzle-plate spacing equation 5.1.1 is solved 
for y and the velocity on the free streamline evaluated by 
equation 5.1.2. Equation 5.1.3 then yields the distribu- 


tion of velocity on the impingement surface. The distribu- 


tions of the mainstream velocity for H = ~, 3.0, 1.0, 0.5 
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Figure 5.1 Two-dimensional jet impinging on a flat surface 
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and 0.25 are shown in Figure 5.2. The curves for H = © and 
H = 3.0 are essentially the same. 

The solution for H = ~ can also be found by consider- 
ing the normal impact of two equal uniform jets, one origin- 
ating at +m and the other at -~. The analysis is described 
in [1] in some detail. When H = © the solution of equa- 
tion 5.1.1 yields,y = "1 and the velocity on the free sur- 
face is Ur = 1 by equation 5.1.2. When y = 1, equation 
5.1.3 reduces toithe velocity distribution described in [1]. 
Also given in [1] is an equation specifying the position of 
the free streamline for H = ~, The free surface is described 


by the equation 
rn} cote fF ex nf 5.1.4 


and is shown in Figure 5.3. 

As shown in Figure 5.2 the rate of change of main- 
stream velocity decreases continually with increasing dist- 
ance from the stagnation point for an infinite nozzle-plate 
spacing. However, as the nozzle-plate spacing is reduced 
the velocity rises more rapidly than linear at a point near 
the lip of the nozzle. This effect may be better seen by 
differentiating the velocity distribution to obtain the 
velocity gradient. Explicit differentiation of equation 
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Figure 5.3 Position of the free surface , two-dimensional jet . 
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Figure 5.4 shows the mainstream velocity gradient distribu- 
tion for nozzle-plate spacings of H ‘=~ 3 3.055 1:0; 0.5 and 
0.25. As H decreases from infinity a peak in the velocity 
gradient distribution curve occurs at a distance from the 
stagnation point equal to the slot half width or at a point 
near the lip of the nozzle. 

The explanation of this phenomena follows the dis- 
cussion in [23]. When the nozzle is placed far away from 
the impingement surface the jet is deflected gradually far 
in advance of the impingement surface. As the nozzle-plate 
spacing is reduced the velocity at the lip of the nozzle 
increases, that is U, > |, and by continuity the velocity 
closer to the stagnation streamline must decrease. Also, 
the flow must turn more sharply in a region which is con- 
tinually decreasing in area, as the nozzle-plate spacing is 
reduced. With decreasing nozzle-plate spacing the stream- 
lines are packed closer together. Because the velocity is 


defined as the derivative of the stream function, the velo- 
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Figure 5.4 Velocity gradient distributions , two-dimensional | 
jet impinging on a flat surface 
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city will be greater if the streamlines are closer together. 
Consequently, the mainstream velocity rises faster than the 
case when H = », 

The finite element program was run for the nozzle- 
plate spacing iota 41350 raha ¢ cOiedtand (0025aertThesacase; 
H = 3.0 was run with both the potential function formulation 
and the stream function formulation. All the other cases 
were run with the potential function formulation only. As 
shown in Figure 5.2 the finite element program approximates 
the distribution of mainstream velocity correctly for the 
various nozzle-plate spacings. The position of the free 
surface, which is part of the solution, is also found accu- 
rately by the finite element method as shown in Figure 5.3. 
Limits of accuracy depend on the number of nodes chosen for 
the approximation and the spacings of these nodes. In 
general, more nodes are required when the velocity is chang- 
ing rapidly or when the shapes of the boundaries are chang- 
ing. capidlytbtafor the ‘gradient of velocity, a cubic spline 
interpolation as described in Appendix 2 is used to numeric- 
ally differentiate the computed velocity distribution. As 
with any numerical differentiation, errors become magni- 
fied. For the case H = 3.0 the finite element method 
accurately approximates the gradient of mainstream velocity 
as shown in Figure 5.4. The velocity varies smoothly and 
numerical differentiation introduces very little error. 
-But for the case of H = 0.25 the velocity gradient has a 


sharp peak and to be able to predict accurately the exact 
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shape of this curve many more nodes than were used for the 


case of H = 3.0 would be required. 


5.2 Heat  lransfer 

Miyazaki and Silberman computed the heat transfer 
near the stagnation point by a finite difference technique. 
The resulting heat transfer curves showing the parameter 
Nup/vRen versus distance from the stagnation point for a 
Prandtl number of ten and nozzle-plate spacings of H = », 
1.5, 1.0 and 0.5 are shown in Figure 5.5. The parameter 
Nup/YRep is related to the temperature gradient at the wall. 


Forming an energy balance at the interface 
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4.1.6 one obtains 


y=0 
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Rearranging and forming a Nusselt number 
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Figure 5.5 Heat transfer near the stagnation point of a two- 
dimensional jet impinging on a flat surface 
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As the nozzle-plate spacing is reduced a peak in the 
heat transfer curve appears near the lip of the nozzle. 
From the discussion of the gradient of velocity in the last 
section it would seem that the heat transfer near the stag- 
nation point is directly related to the gradient of the 
mainstream velocity. To show this relationship the simple 
case of Hiemenz stagnation point flow discussed in White 
[18] will be used. 

The stagnation point flow of an infinite fluid or 
Hiemenz flow is one of the few flows having an exact solu- 
tion of the complete Navier-Stokes equations. The solution 
is described in [18] and involves a similarity parameter 


n, defined as 


: n= eye 
and a non-dimensional form of the stream function F(n) 
defined by 
dU. io8 
yp = Kora V Ga) 5.2. 
dx 0 


where Gan is the mainstream velocity gradient at the 


stagnation point and v is the kinematic viscosity. The 


mainstream velocity for this stagnation point flow rises 
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linearly with the distance from the stagnation point or 


Us Bx Sere 


where the parameter 8 may be thought of as the velocity 


gradient at the stagnation point. 
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Substitution of equations 5.2.3 and 5.2.4 into the Navier- 
Stokes momentum equation yields an ordinary differential 
equation for F. 
3 2 
ee FS T= 0 BeeeT 
dn dn a 
with boundary conditions 
F(0) = F (0) = 0; F'() = 1 5.2.8 


the solution,0¢.5..247 .with boundary GOnd1 Lions. oo2.o 1S 
available in either [17] or [18] in tabular or graphical 


form. 


Now consider the energy equation and a dimensionless 
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where T, is the ambient fluid temperature and To the surface 
temperature. The energy equation, neglecting viscous dis- 


Sipation becomes 


c 
dO do 
oe a 


using equations 5.2.3, 5.2.4 and 5.2.9. The boundary condi- 


tions are 


e(0) = 0; Q@(o) = ] Bae ed 


Equation 5.2.10 with boundary conditions 5.2.11 can be 


directly integrated to give 


n n 
{exp |- Prinietokds 7 di} 
a(n) = 2 - 5.2.12 
{exp |- Pr f Fds| dn 
0 0 


The heat transferred at the wall is then given by 
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Forming a Nusselt number results in 


eshdD . 
Nup i ae Bec ed 
where G(Pr) is given by 
foe) n o7 
G(Pr) = ¢f exp(- Pr f  Fds)dn ier 15 
0 0 


For Prandtl number equal to ten G(Pr) is 1.3388 from [18]. 
Using the non-dimensionalizations of equation 2.1.12, equa- 


tion 5.2.6 becomes 


dU dU 
Pre RO Si 


and evaluating equation 5.1.5 at the stagnation point, 


U. = 0, gives 
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5.2.14 and using the value of G(Pr) for a Prandtl number 
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The solutions of equations 5.1.1, 5.1.2 and 5.2.18 for 
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Seeake 


various nozzle-plate spacings and a Prandtl number of ten 


are given in Table l. 


Table 1. 
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.00000594 


000263 
. 004080 
. 044880 
-474190 
«451139 
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.000771 
.016355 
.065911 
-234042 
.902776 
423336 


The values of Nupy/vRep for H = = 


given in Table 1 are shown at the stagnation point on 


Figure 5.5. 


Exact solutions of a two-dimensional 
jet impinging on a flat surface. 
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by the finite difference and integral boundary layer methods 
at the stagnation point. Therefore, Heimenz flow is valid 
at the stagnation point for a finite jet. Because the velo- 
city in Heimenz flow is linear by equation 5.2.5, the 
gradient of velocity will be constant. By equation 5.2.14 
this implies a constant Nusselt number and constant heat 
transfer in the region of the stagnation point. As the 
Nusselt number varies as the square root of the velocity 
gradient, by equation 5.2.14, the Nusselt number and heat 
transfer are not constant for the case of the finite jet. 

As seen by this analysis the peaks in the gradient of velo- 
city distribution are directly related to the peaks observed 
in, the. heat, transfer curves« 

To compare the integral boundary layer technique of 
Chapter IV with the finite difference method of [23], the 
integral boundary layer solution was computed for nozzle- 
pi ahex~spacings, ofate te 1,5.) beOse0e5eand: 0.25. for the case 
of no melting and a Prandtl number of ten. The non- 
dimensional heat transfer Nup/YRep is modeled correctly 
by the integral boundary layer solution for nozzle-plate 
spacings of 1.0 and 0.5 as shown in Figure 5.5. The seem- 
ingly oscillatory nature of the integral boundary layer 
solution is due to the inaccuracies in calculating the velo- 
city gradient by cubic spline interpolation. 

The discrepancy for the case of H = 1.5 is explained 
-as follows. Miyazaki and Silberman [23] concluded that the 


velocity distribution curves for H = ~ and H = 1.5 were 
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identical. As shown by equation 5.2.14 the Nusselt number 
varies as the square root of the streamwise velocity 
gradient at the stagnation point. Evaluating equation 


She. bhufonrh =sceandchia=ihad result «in 


dU, 1/2 
(a5) = 0.886 for H = = 

5.2.19 
pase note be fonatieies 
dx ’0 


The curve presented by Miyazaki and Silberman [23] 
for H = ~, 1.5 probably used the velocity distribution of 
equation 5.1.3 for H = ~. The finite element solution used 
a nozzle-plate spacing of H = 1.5. Therefore the finite 
element solution should be 0.914/0.886 times greater than 
the solution presented in [23]. This is actually the case 
as :sfrown “inw"Figure’ 57%. 

The solution for'a nozzle-plate spacing of H = 0.25 
was attempted but because of the large velocity gradients 
present, as shown in Figure 5.4, the boundary layer solution 
failed. As was pointed out in Chapter IV the range of 
pressure gradient parameter for which the solution is valid 
foeetces Ns she oes Ded The lower limit corresponds to 
separation and the upper limit to the case where the velo- 
city profiles pop [18]. When the nozzle-plates spacing is 
0.25 the large velocity gradients imply that a large pres- 


sure gradient is present. 
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In both cases, for too small a pressure gradient 
A < -12, and for too large a pressure gradient hey p> he 
(1 + 2r9) the integral boundary layer solution is not valid. 
The boundary layer equations 4.1.2 and 4.1.3 no longer apply 
and to solve the problem in the regions A < -12 or A > 12 
(1 + 29) another method of solution must be considered. 

The errors introduced by the finite element method 
and integral boundary layer solution are caused by many 
factors. The number of nodes chosen affects the accuracy 
of the calculated velocities. In general, the more nodes 
used to approximate the flow field, the greater the accur- 
acy will be but more computer time is required. Therefore 
there is a trade off between accuracy and computer time used. 
From comparison of the finite element generated velocity 
distribution and the exact solution equation 5.1.3 the 
approximate accuracy for all cases of nozzle-plate spacings 
is less than +5 percent. The same number of nodes were 
used in each case. The greatest errors occur when the velo- 
city is changing very rapidly, that is when H = 0.25. The 
accuracy of the velocity gradient calculation depends prim- 
arily on the accuracy of the cubic spline numerical differ- 
entiation procedure. Again comparing the finite element 
solution to the exact solution of equation 5.1.5 in Figure 
5.4 the approximate accuracy is +20 percent. Again, greater 
accuracy is achieved when the velocity gradient does not 
change very rapidly, that is for H = 3.0. 


In the calculation of the heat transfer there are 
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two main sources of error. One is the error introduced by 
the approximation of integrating the boundary layer equa- 
tions and using a fourth order polynomial for the velocity 
and temperature profiles. The other error is introduced 

by the approximation of the velocity gradient in the cubic 
spline approximation. With the Nusselt number being pro- 
portional to the square root of the velocity gradient and 
an error in the velocity gradient of +20 percent the maximum 
error in the Nusselt .number is approximately. +10.6 percent. 
Because of the averaging quality of the integral solution 
an approximate guess to the error in the calculation of the 


Nusselt number is between 10 and 15 percent. 
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CHAPTER VI 


AXISYMMETRIC JET IMPINGEMENT 


As mentioned previously there are two distinct cases 
of jet flows. Both cases apply to axisymmetric and two- 
dimensional jets. The first occurs when a gaseous jet 
issues into a gaseous medium or a liquid jet issues into 
a liquid medium. Here, entrainment of the ambient fluid 
by the jet greatly affects the velocity profile of the jet. 
The second case occurs when a liquid jet issues into a 
gaseous medium. This jet is virtually unaffected by the 
ambient fluid and a free surface forms between the liquid 
jet and the surroundings. 

The impingement characteristics of these two jet flows 
are slightly different. The flow regions for the case of 
the entrained jet are well known. Near the nozzle exit a 
potential core region exists, surrounded by a mixing region 
between the jet and the ambient fluid. A few jet widths 
or diameters downstream the mixing region has spread inwards 
to engulf the potential core. Beyond this point the mixing 
continues and the velocity distribution adjusts itself so 
as to conserve linear momentum. In this region the velo- 
city profiles are self-similar. When the jet impinges on 
a surface, in the stagnation region, a strong essentially 


inviscid jet-surface interaction produces a change in flow 
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direction [24]. Far from the stagnation point a wall jet 
forms and the velocity profiles are again self-similar. 
A transition zone exists between each of the above regions. 

The flow regions for the case of the free surface jet 

are slightly different. The velocity profile present at the 
exit of the nozzle persists until the stagnation region 
begins to change the kinetic energy of the jet into pres- 
sure energy. In the impingement region and resulting flow 
over the surface Watson [25] recognizes four distinct flow 
regions. These are shown in Figure 6.1 and are described 
below. 

1. A stagnation region, where the velocity on the 
surface rises rapidly from zero at the stagnation 
point to the free streamline velocity at a dist- 
ance from the stagnation point equal to approxi- 
mately the nozzle diameter. When x' << D the 
flow can be described by the stagnation flow of 
an infinite fluid where the mainstream velocity 
rises linearly with x. This is the Hiemenz 
stagnation point flow described in [17] or [18]. 

2. \A flat plate region at a greater distance x' where 
there is no pressure gradient and the flow is 
similar to Blasius flow over a flat plate. 

3. A transition region where the boundary layer 
includes the whole flow field and the velocity 
profile changes from the flat plate profile to 


the fully developed similarity profile. 
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VELOCITY PROFILE 
AT NOZZLE EXIT 
UNCHANGED 


IMPINGEMENT SURFACE 
ALTERS VELOCITY 
PROFILE 


TRANSITION 
REGION 


STAGNATION FLAT PLATE SIMILARITY 
REGION REGION REGION 


Figure 6.1 Flow regions of an impinging jet 
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4. A similarity region where the fully developed 
Similar profile persists. A hydraulic jump will 
eventually occur and the thickness of the liquid 
on the plate suddenly increases. 

The above regions apply only to laminar flow. As 
region 1 has a strong favorable pressure gradient, transi- 
tion to turbulence is not likely to occur there. A more 
likely place for transition is in region 2 where the pressure 
gradient is essentially zero. Of course transition to 
turbulence will depend on a number of factors, jet Renolds 
number, free stream turbulence, and length scales and sur- 
face roughness. Only laminar boundary layers are considered 


in this study. 


6.1 Previous Work 

The main problem in determining the flow field char- 
acteristics of axisymmetric jets is that the stream function 
and the potential function both do not satisfy Laplace's 
equation and thus the methods of complex variables are not 
useful. Therefore, numerical and experimental methods must 
be attempted. Leclerc [26] used an electrical analogy to 
solve for the potential flow of an axisymmetric jet imping- 
ing on a flat plate. Brady and Ludwig [27] computed the 
velocity on a flat plate due to an axisymmetric jet by 
experimentally measuring the pressure distribution on the 
surface. Their interest was in the ground flow characteris- 


tics of VTOL and STOL aircraft jets. Scholtz and Trass 
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[28] used a separation of variables technique and a trun- 
cated series solution to solve for the flow field in a non- 
uniform impinging jet. Sarpkaya and Hiriart [16] used the 
finite element method of [11] and also Belotserkovsky's 
integral method to solve for the potential flow associated 
with aerodynamic jet thrust reversers. 

Other investigators [29]-[39] have made experimental 
and theoretical investigations into the flow field and heat 
transfer of impinging axisymmetric jets. The effects of 
entrainment are important and some of the parameters and 
effects dealt with are nozzle-plate spacing, turbulence 
levels and length scales, size of jet and size of target 
surface, jet Reynolds number and Prandtl number. The pres- 
sure distribution measured on the surface is important for 
determining the velocity on the impingement surface. 

In this chapter the flow and heat transfer characteris- 
tics of an axisymmetric laminar jet impinging on a surface 
will be investigated. Two effects were modelled by the 
finite element method and integral boundary layer analysis. 
These are the effect of nozzle-plate spacing and the effect 
of surface curvature. Also the shear stress on the imping- 
ment surface is determined. The effects of melting on the 


heat transfer will be considered in the next chapter. 


6.2 Flow Field of Axisymmetric Jet Impingement 
The finite element program was run for an axisymmetric 


jet impinging on a flat surface for nozzle-plate spacings 
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Open) = 10225) Oo eve anaes.0. The radial surface velocity 
distributions are shown in Figure 6.2. The velocity distri- 
butions for the cases of 1.0 and 3.0 are essentially the 
same and therefore the case of H = 1.0 corresponds to the 
case of infinite nozzle-plate spacing. This is contrasted 
with the two-dimensional case where the infinite nozzle- 
Bate spacing wis modelled by H = 3.0. —In the region x < 0.25 
the velocities for all three nozzle plate spacings are 
approximately equal, although the velocity gradients are 
very different. This is shown in Figure 6.3 where the velo- 
city gradients for H = 3.0, 1.0, 0.5 and 0.25 were calcu- 
lated by the finite element method and cubic spline inter- 
polation technique. As was present in the two-dimensional 
case errors occur due to numerical differentiation of the 
computed velocities. Peaks occur in the gradient of velo- 
city curves at a radial distance near the lip of the nozzle. 
Similar to the two-dimensional case, these peaks are related 
to peaks which occur in the heat transfer curves. A differ- 
ence from the two-dimensional case is that a peak occurs 

in the i¢ase,H sa1.0, © for the axisymmetric jet. There is 
no peakwin®the gradient of velocity cumve for H = 93.0, = 

in the two-dimensional case. The peak in the axisymmetric 
case is due to the rapidly decreasing flow area near the 

lip of the nozzle so that by continuity requirements, the 
velocity must increase more rapidly than in the two-dimen- 
sional case. 


The velocity distribution predicted by Leclerc's 
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FINITE ELEMENT 
© H=3.0,1.0 


O H=0.5 


m H=0.25__ 


0.0 0.2 0.4 0.6 0.8 1-0 1.2 


X 


Figure 6.3 Velocity gradient distributions , axisymmetric 
jet impinging on a flat surface 
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electrical analogy [26] for H = ~ and by the experimental 
method of Brady and Ludwig [27] for H = ~, 0.5 and 0.25 

are also shown in Figure 6.2. Agreement between the finite 
element method and the experimental methods is quite good 
for the case of H = 1.0 and H = 0.5, but the accuracy 
suffers when H = 0.25. All three cases had the same number 
of nodes in the finite element mesh and the decrease in 
accuracy is due to the increase in velocity gradient. A 
large change in velocity requires more nodes for accurate 
representation. 

_The same argument for the increase in velocity grad- 
ient with a decrease in nozzle-plate spacing as given for 
the two-dimensional case also applies for axisymmetric jet 
impingement. 

To consider the effect of a curved impingement sur- 
face the finite element program was run with a nozzle-plate 
spacing of H = 3.0 for an axisymmetric jet impinging on 
curved paraboloid surfaces. The paraboloid cross-sections 


were of the form 


for various values of the parameter a, (a = 0, 0.05, 0.1, 
0.2, 0.5 and 1.0). The coordinate directions w and r along 
with the resulting velocity profiles are shown in Figure 
6.4. The radius of curvature R and curvature K are found 


by the formula 
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The smallest radius of curvature and largest curvature occur 
at r = 0, the stagnation point. Table 2 shows the radii 
of curvature and curvature for the family of paraboloids 


considered. 


Table 2. Radii of curvature and curva- 
ture for a family of parabo- 
loids given by w = ar¢ at r = 0. 


a K R 
0 0 co 
0.05 0.4 10 
0.1 0;2 5 
0.2 0.4 236 
O25 1e0 1.0 
1.0 250 0.5 
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For the boundary layer equations to be valid the 


relationship ae << 1 must hold [18]. Using the non-dimen- 


Sional variables of equation 4.1.6 and the relationship 


R' Sy “1/2 


R = an the curvature criteria becomes “Re (Rep eee les 


If this criteria is not fulfilled then the complete Navier- 
Stokes equations including curvature effects must be con- 


§ 
sidered [45]. 1f (GUMS << 1 is valid, then the effects 


of a change in surface shape are felt through the external 
pressure gradient the surface imposes on the flow. For flow 
on the concave side of the paraboloid the pressure gradient 
is favourable and tends to delay separation. For flow on 
the convex side, the pressure gradient is unfavorable and 
tends to increase the chances of separation occurring. Only 


the flow on the concave side of the paraboloid was considered 


6 
assuming that the curvature criteria, — Rene os << 1, was 


satisfied. 


6.3 Heat Transfer of Axisymmetric Jet Impingement 

With the velocity distributions of the previous sec- 
tion the boundary layer solution described in Chapter IV 
may be used to find the heat transfer characteristics of 
the impinging axisymmetric jet. As outlined in the previous 
chapter concerning the two-dimensional jet the Nusselt 
number was found to be proportional to the square root of 
the velocity gradient at the wall. The analysis for the 


axisymmetric jet is the same as both the two-dimensional and 
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axisymmetric stagnation point flows are of the same type. 

If boundary layer equations are considered, they both belong 
to the class of Faulkner-Skan wedge flows [18]. The diffi- 
culty for axisymmetric jet impingement lies in finding the 


parameter B in the expression 
28 
U_ = a x 6.3.1 


As the exact solution is available for the two-dimen- 
sional jet, the parameter 8 was found as shown by equations 
Dc and 524s. 

The governing equation for axisymmetric Heimenz stagna- 


tion point flow is [18] 
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similar to that for the two-dimensional case. The only 
difference is the presence of the factor 1/2 in the above 
equations. The boundary conditions are the same as for the 
two-dimensional case, that is F(0) = F'(0) = 0 and F'(~) = 
1. The energy equation and thermal boundary conditions are 
identical with those of the two-dimensional case, that is 
equations 5.2.10 and 5.2.11. The solution is again equation 
5.2.12 but the non-dimensional stream function F is obtained 
from the solution of equation 6.3.2. The expression for 

the Nusselt number is given by equation 5.2.14 and is 


repeated here for reference. 


ao pele 
Nup = DiS G(Pr) Oro 45 


The value of G(Pr) for axisymmetric stagnation point flow 
and a Prandtl number of 13.7 is 1.3296 as given by White 
[18]. This compares with the two-dimensional value of 
1.4557 for a Prandtl number of 13.7. The physical proper- 
ties, assumed constant for this problem, must be calculated 
at some reference temperature. Three possible choices are 
the wall temperature, the film temperature and the jet bulk 
temperature. The wall temperature, O0°C, was chosen as the 
reference temperature for calculating the physical proper- 
ties. A discussion on the applicability of this assumption 
is given in the next chapter. Using the non-dimensionaliza- 


tion of equation 2.1.12, equation 6.3.1 becomes 
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similar to equation 5.2.16. Substituting equation 6.3.6 


into equation 6.3.5 and rearranging results in 


Nup elk oes 
= (A repre) WNL Nu ea 6.3.7 
/Re, ita 


As no exact solution for axisymmetric jet impingement 
exists the values of 8, the velocity gradient at the wall, 
may be obtained from the finite element and cubic spline 
interpolation. The values of the velocity gradient at the 
stagnation point from Figure 6.3 and the heat transfer para- 
meter from equation 6.3.7 are shown in Table 3 for various 


nozzle-plate spacings. 


Table 3. Axisymmetric jet impingement heat 
transfer on a flat surface. 
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The values of Nup/vRep from Table 3 at the stagnation 
point for’ H = ©, 1:0, and 0.5 are shown in Figure 6.5. 
Also shown are the heat transfer curves generated by the 
boundary layer analysis. The case of H = 0.25 was attempted 
but the velocity gradients were too large. Large velocity 
gradients imply a large pressure gradient parameter A, and 
the solution is not valid for large A as previously des- 
cribed. The scatter in the data is due to the numerical 
differentiation procedure of the cubic spline interpolation 
approximation coupled with the integration procedure. Peaks 
in these heat transfer curves correspond to the peaks in 
the velocity gradient distribution curves of Figure 6.3. 

The trends in heat transfer on curved surfaces can 
be inferred from the velocity distribution curves in Figure 
6.4. The heat transfer for the curve w = 0 is same as 
shown in Figure 6.5 for H = 3.0, 1.0. As the parameter 
a increases, the velocity gradient at the stagnation point 
decreases as shown in Figure 6.4. Therefore, the stagnation 
point heat transfer will also decrease with increasing 
curvature. This trend is the same as concluded by Cheng 
and With@ane boa). ) For thercases,.a°=.0.05.) 0-1) amte0, 22 the 
heat transfer curves will only slightly differ from the 
flat surface case a = 0. This is evident because the velo- 
city distributions for these four cases are similar. For 
the cases a = 0.5, 1.0 the velocity gradient in the region 
near the stagnation point is very small. Heat transfer by 


conduction will predominate over convective heat transfer 
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0 H=3.0 ,1.0 
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Figure 6.5 Heat transfer near the stagnation point of an 
axisymmetric jet impinging on a flat surface 
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in the region of the stagnation point for these two cases. 
Also of interest is the shear stress generated on 

the surface due to the impinging jet. The shear stress 

will be important if erosion effects are present in the 

transfer mechanisms. A possible situation where this would 

happen is the impingement cutting of a frozen sand or gravel. 


Define a shear stress coefficient as 


i du’ 
i. ly ele 
x = D 6.3.48 


oV 


Introducing the non-dimensional variables of equation 4.1.6 


results in 


Equation 4.4.1 is used to evaluate the velocity gradient 

in the direction normal to the wall. Figure 6.6 shows the 
shear stress distribution on the impingement surface as 
calculated by the integral boundary layer analysis. From 
equation 4.4.1 and equation 6.3.9 the shear stress parameter 
is proportion to the ratio U/S |. The shear stress para- 
meter is zero at the stagnation point as the free stream 
velocity is zero there also. For the region near the 

the boundary layer thickness, is con- 
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Figure 6.6 Shear stress on a flat surface due to axisymmetric 
jet impingement 
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stant and U rises from zero as shown in Figure 6.2. There- 
fore, the shear stress parameter rises from zero in a 
manner similar to the rise in the velocity. The smaller 
nozzle-plate spacing produces a larger shear stress para- 
meter because a greater rise in free stream velocity is 
experienced for smaller nozzle-plate spacings as shown in 
Figure 6.2. As the radial distance increases the main- 
stream velocity becomes constant and the boundary layer 
thickness increases. These two effects combine together to 
produce the peaks in the shear stress parameter curves near 
the lip of the nozzle. The subsequent decrease in shear 
stress is due to the continual increase in boundary layer 
thickness while the mainstream velocity remains constant. 
The combined effect of shear stress and melting will 


be discussed in the next chapter. 
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CHAPTER VII 


MELTING AND EROSION CAUSED BY AXISYMMETRIC 
JET IMPINGEMENT 


The melting and erosion of a frozen material by the 
action of a water jet presents a very complex problem. 
Because of this complexity, investigations on melting by 
impinging jets have been experimental in nature. In this 
chapter the melting and erosion of frozen material is 
analyzed by the methods of the previous chapters. The 
effect of melting on the heat transfer of an axisymmetric 
jet impinging on an ice surface is shown. An attempt to 
model the melting of a block of ice results in the steady 
state heat transfer rate as a function of Stefan number. 
To obtain a physical perspective of the melting process, 
photographs were taken at various time intervals during the 
medting of a Diock Of ice.. Finally the effect of shear 
stress and erosion on the cutting action of the jet is 


discussed. 


7.1 Previous Work 

Recently experimental investigations on the melting 
of ice with a water jet have been made. Savino, Zumdieck 
and Siegal [40] have studied the effect of freezing and 


melting on the heat transfer coefficient at a stagnation 
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point and found that the phase change had little effect 

on the heat transfer coefficient. However, they only con- 
sidered small temperature differences between the jet and 
the ice and the Stefan number was small in their study, 

Ste < 0.14. Yen [41] experimentally studied the heat trans- 
fer characteristics of a bubble-induced water jet impinging 
on an ice surface. In this case the water jet is submerged 
and entrainment is an important factor. As previously 
mentioned Mellor [42] has presented a number of potential 
applications for the cutting of ice with continuous high 
pressure jets. Two recent experimental studies by Yen and 
Zehnder [43] and Gilpin [44] have shown the effects of melt- 
ing a block of ice with a water jet. In [43] the mass of 
ice removed by the ablation process was the primary mea- 
sured quantity. They found that the weight-loss versus 

time curves were essentially linear implying a constant melt- 
ing rate. Gilpin [44] recognized two distinct regimes of 
melting heat transfer, characterized by the shape of the 
cavity in the ice produced by the impinging water jet. The 
first was a smooth cavity shape occurring at Reynolds 


numbers in the range, 2 x 10? 


<Repn’< 4°x 107. Heat trans- 
ferred during steady state cutting of the ice was found to 


be correlated by 
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The exponent 0.65 indicates a heat transfer rate slightly 
greater than that for purely laminar flow. Smoothness in 
cavity shape implies that a laminar boundary layer is pre- 
sent and the increased rate of heat transfer is due to the 
turbulent free stream of the jet. The second case was a 
rough cavity shape indicating highly turbulent flow. The 
Reynolds number range for this case was .5 x 10° < Rep iS 
ax 703 Obviously either mode of heat transfer can exist 
for a range of Reynolds numbers and the mode which does 
exist depends on the experimental conditions at the time 


the experiment takes place. For the rough case the heat 


transfer was correlated by 
Nu N= OND @uRe ce” Pale 


The exponent 0.8 is characteristic of turbulent forced con- 


vection heat transfer. 


7.2 Melting Ice with a Water Jet 

The boundary layer solution was run for an axisymmetric 
jet impinging on a flat ice surface. The heat transfer 
distribution near the stagnation point for a range of Stefan 
numbers is shown in Figure 7.1. Since the water temperature 


is limited to the range 0°C to 100°C the Stefan number, 
ieee ft 


is limited to the range 0 to 1.25. The nozzle- 
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Figure 7.1 Heat transfer from an axisymmetric jet to a flat 
surface , effect of melting 
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plate spacing used was H = 3.0. For increasing Stefan 
number the heat transfer parameter Nup/vYRep, decreases. 
The shape of the heat transfer curve is independent of the 
Stefan number but shifts downwards with increasing Stefan 
number. The decrease in heat transfer with increasing 
Stefan number may be explained as follows. From equation 
5.2.2, the heat transfer parameter is equal to the non- 
dimensional temperature gradient at the wall. Obviously, 
as the bulk temperature of the jet is increased more ice 
will be melted and a thicker layer of near freezing water 
will be formed. With the thicker layer of colder water, 
the temperature gradient at the wall decreases. A shielding 
of the ice surface by the newly melted water layer takes 
place. 

The finite element method and integral boundary layer 
solution were used to model the continuous melting of a 
block of ice by a water jet. The solution procedure pro- 
ceeded as follows. First, the axisymmetric jet flow field 
on a flat surface was solved by the finite element method. 
The integral boundary layer solution was then solved to 
obtain the distribution of melting velocity, Vo» along the 
ice surface. If the position of the ice surface before the 
melting started was z+ then the new position of the ice 


surface after a non-dimensional time At is 
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using the coordinate system described in Figure 7.2. An 

adjustment of the new ice surface position now takes place 
to ensure that the finite element mesh does not become too 
distorted from it's original shape. If Z+0 is the original 
position of the stagnation point and Z(t+At)0 the new posi- 
tion of the stagnation point after melting has taken place, 
then the melted ice surface is adjusted according to the 


transformation 


Za = Zt+at ~ (2(t+at)o ~ to) oo! 


where Zy gives the adjusted ice surface shape. This trans- 
formation ensures that the stagnation point remains at the 
same vertical distance from the lip of the nozzle. The 
transformation for the first time step is shown schematic- 
ally in Figure 7.2. With this transformation the origin 
of the coordinate system remains at the stagnation point. 
Physically the adjustment of the ice surface may 
correspond to two entirely different situations. As the 
distance between the nozzle and the stagnation point remains 
constant the first situation occurs when the nozzle moves 
with the same velocity as the melting surface. The relative 
velocity between the water jet aia the-ice surface is V, 
the dimensional water jet velocity. The second case occurs 
if the nozzle remains stationary and the ice surface melts 


away from the nozzle. Although the distance between the 
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Figure 7.2 Adjustment of the melted surface 
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nozzle and the stagnation point is physically increasing, 
a nozzle-plate spacing of H = 1.0 models all cases in the 
range 1.0 < H < ~ as previously discussed. The initial 
nozzle-plate spacing for this case was H = 3.0 and the 
adjustment of the ice surface does not alter the potential 
flow analysis. However, the relative velocity between the 
water jet and the ice surface is now V - Vo: The program 
was run assuming the relative velocity was V, that is, the 
first case was considered. As will be shown later the 
magnitude of the melting velocity is so small that the 
distinction between these two cases is negligible. 

The complete computer program was run for Stefan 
numbers of Ste = 0.2, °0.33, O772 Mand 1217. A typicak time 
history profile of the ice interface is shown in Figure 
7.3 for a Stefan- number of 0.33.e@,A time step of @t = 1 
was used for this run. The non-dimensional times for the 
ice interface profiles shown are t = 0, 5, 10, 20 and 30. 
The ice profile is atnost flatefor the region 0 <x < 0.75. 
As time increases the ice profile becomes steeper in the 
region 1.0 < X < 1.75. The program was stopped after 30 
iterations because separation was predicted. Separation 
occurs when the pressure gradient parameter A becomes less 
than -12 as previously discussed. 

The assumption of quasi-steady state flow made in 


Chapter IV may now be verified. Consider a typical flow 
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The Reynolds number for this flow is approximately 8.6 x toe: 
For a Stefan number of 0.33, the non-dimensional melting 
velocity is approximately ees 0.04. Therefore the actual 


melting velocity for this problem is 


or only about 0.04 percent of the incoming jet velocity. 
Because the melting velocity is such a small percentage of 
the incoming jet velocity the flow may be assumed to be at 
steady-state. 

Also because of this small percentage the difference 
between the two cases of relative velocity will be very 
small. Therefore, the jet impingement problem may be 
modelled by the case where the nozzle moves with the same 
velocity as the melting ice surface. 

As shown in Figure 7.3 the shape of the ice profile 
in the flat region is constant for the last ten iterations. 
The constant shape implies a constant melting velocity along 
the whole flat region of the ice profile. Steady state 
melting was assumed when the ice profile in the flat region 
had reached a constant shape. The melting velocity Vos is 
then used to calculate the steady state heat transfer rate. 


From equations 4.1.11 and 5.2.2 the heat transfer 


parameter is given by 
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For the case of no melting both the melting velocity 
and the Stefan number are zero and equation 7.2.2 is indeter- 
minate. Figure 7.4 shows the steady state heat transfer 
parameter versus Stefan number. The region near Stefan 
number equal to zero is shown by a dotted line to indicate 
that the steady state results for this region are not known 
because of the indeterminateness of equation 7.2.2. Also 
shown are the stagnation point values from Figure 7.1. As 
the Stefan number increases, the difference between the 
initial stagnation point melting velocity and the steady 
state melting velocity increases. This trend leads to the 
assumption that the steady state heat transfer rate is the 
stagnation point heat transfer rate for a Stefan number 
equal to zero. Therefore, for small Stefan numbers the 
initial stagnation point melting velocity may be used to 
calculate steady state heat transfer rates but at larger 
Stefan numbers differences of up to about ten percent would 
result if stagnation point values were used instead of the 
steady state value. 

The next problem was that of deciding on a reference 
temperature. As pointed out in the derivation of the 
boundary layer equation, the physical properties were 
assumed to be constant. The water temperature varies from 


the melting temperature at the wall to the bulk temperature 
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Figure 7.4 Steady-state heat transfer due to axisymmetric 
jet impingement 
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of the jet in the free stream across a thin boundary layer 
and the problem of which temperature in this range is to 

be used for calculating the physical properties arises. 

This is a problem as the viscosity of water varies greatly 
over the range of Stefan numbers considered [46]. The heat 
transfer curves for references temperatures based on the 
wall temperature, the film temperature and the bulk tempera- 
ture are shown in Figure 7.5. A significant variation in 
the heat transfer parameter is noticed depending on which 
reference temperature is used. To determine which reference 
temperature is best suited to describe this heat transfer 
process some experiments should be made and the heat trans- 
fer calculated based on each of the three reference tempera- 
tures. A comparison could then be made with Figure 7.5. 
Another possible solution would be to attempt a variable 
property numerical solution for this jet impingement heat 
transfer problem. 

To gain some physical insight into this forced convec- 
tion melting heat transfer problem a series of simple experi- 
ments were conducted. The experimental apparatus is shown 
schematically in Figure 7.6. A jet of water was vertically 
directed at a block of ice for varying amounts of time. 

The water was pumped from a constant temperature tank and 
the temperature of water and flow rate were recorded. After 
completion of run, the block of ice was removed from the 
-apparatus and the block was cut in half to reveal the pro- 


file of the hole produced by the impinging water jet. A 
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Figure 7.5 Effect of reference temperature for calculating. 
physical properties 
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Figure 7.6 Experimental apparatus 


me Oe 


ied au 
i 
~ =) 
2) 
“Ss 
‘ 
a x 7 


ri m 


asmienutT Mua 
asTaMWOla 


oe ~~ 
39 


a 


- toe yy ai i! 


> 
a j 
at es} eee | pon 


ae 


wor 
- 
= 


ee 
— 
ea 
Cee a= ae 
a 


- 


erste 


138 


photograph of the resulting ice profile was then taken. 

The time history of the ice interface was obtained by allow- 
ing the jet to impinge on the ice surface for various times. 
A 1.27 cm diameter jet of water at a temperature of 25°C, 
Stefan number of 0.31, was directed at an initially flat 
block of 1cemGOmestshmnesenG fOr Ose t Once dOwandeSOmseconds. 

The resulting photographs are shown in Plate 1. The jet 
velocity was approximately 0.76 meters/second. Using the 
wall temperature as reference temperature this corresponds 
to a jet Reynolds number of Rep = 5.4 x 10”. To compare 


the photographs with the curves in Figure 7.3 a reference 


time must be established. Knowing 
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For the present case the reference time is 1.22 seconds. 
The times considered in the experiment then correspond to 


non-dimensional times of 4.9, 8.2, l6s44tems and 65.6. A 
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1c) t' =20 seconds, t = 16.4 


1d) t’ =30 seconds, t = 32.8 


te) t' = 80 seconds, t = 65.6 
Plate 1. Time history of a block of ice melted 


by a water jet. D = 1.27cm, V = 0.76m/sec 
Re, = 5.4X103, Ste = 0.31. 
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discussion of the photographs follows. 

Plate la. A very shallow hole is present. The shape 
of the ice interface looks like that in Figure 7.3 for non- 
dimensional time of 5 in the region X < 1.5. Beyond this 
point, transition to turbulent flow is present, increasing 
the heat transfer and changing the shape of the ice inter- 
face. 

Plate 1b. A deeper hole is shown. Note the flatness 
of the region near the stagnation point and the sudden rise 
in profile shape. The region of turbulent flow can be seen 
very clearly. 

Plate 1c. An even more sharply defined laminar type 
region. A flat region still exists near the stagnation 
point. 

Plate Id. A very sharp distinction at the end of the 
laminar flow region. The flow is separated and continual 
application of the jet just deepens the cavity and increases 
the steepness of the cavity wall. The stagnation region is 
still very flat. 

Plate le. The flow has melted away the sharp lip 
shown in Plate 1d and has reattached to the block of ice. 
The formation of another lip, and then separation is pre- 
dicted to occur. This pattern of separation and reattach- 
ment is predicted until the hole has formed nearly vertical 
walls and the flow is continually separated. 

From comparison between the photographs and the pre- 


dicted results of Figure 7.3 the conclusions are: 
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1. That a flat region exists near the stagnation 
point. 

2. That the computer program is valid only in the 
region 0 < X < 1.5 because transition to turbu- 
lence is likely to occur beyond this point. 

3. That the computer program predicts separation 
accurately. Also the computer analysis is not 
valid after separation has been predicted to occur. 

The effect of shear stress as well as heat transfer 

may be important in some applications where a jet is used 
to cut a frozen medium such as frozen gravel or sand. 

Assume now that in addition to the heat transfer the removal 
rate of the material is directly proportional to the shear 


stress so 


i = A aE 
0 pL oO 


Toa 4 
Yay 


i} 
Ss 


where c is a proportionality factor. Introducing the non- 
dimensionalizations of equations 4.1.6 and equation 6.3.8 


results in 
2 oles oF ec a 


-A problem arises because the value of the factor c¢ is unknown. 


To determine the effect of shear stress the program was run 
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for the case i (Rey =! and the resulting shape of the inter- 
face after one melting time step is shown in Figure 7.7. 
The contribution of the shape due to melting is shown for 
comparison. The shape of the interface in Figure 7.7 
resembles the inverse of the shear stress coefficient curve 
of Figure 6.6 which was expected considering the form of 
equation 7.2.5. 

| The exact value of the parameter c is unknown but the 
effects of varying the parameter 5 YREp are possible to 
predict with the present program. It is recommended that 
an experimental program be set up to determine what factors 
affect the parameter c and what it's value might be for 


certain materials. 
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CHAPTER VIII 


CONCLUSIONS AND RECOMMENDATIONS 


This thesis was motivated by a desire to be able to 
predict the melting characteristics of a water jet imping- 
ing on an ice surface. With this in mind an analytical 
solution for laminar jet impingement has been derived. 

Both two-dimensional and axisymmetric flows are dealt with. 
The solution is divided into two parts, a potential flow 
problem which is solved by the finite element method and a 
boundary layer problem, solved by the Karman-Pohlhausen 
integral method. 

the. potential flowssoimtion 1s.\valid tor all jet 
impingement free surface type flows where gravity effects 
can be neglected. The physical effects dealt with are the 
distance from the nozzle-exit to the impingement surface 
and the impingement surface shape. From the potential flow 
solution the position of the free surface and the velocity 
distribution on the impingement surface used for the 
boundary layer analysis, are found. 

The boundary layer solution is valid for laminar 
boundary layers on arbitrarily shaped surfaces. However, 
if the curvature of the surface becomes too large the bound- 
ary layer equations are no longer applicable and the com- 


plete Navier Stokes equations, including curvature effects, 
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must be considered. Another restriction on the applicabi- 
lity of the solution is imposed by the integral solution 
technique. The Karman-Pohlhausen integral technique is 
limited by extreme values in pressure gradient which may be 
present in the flow. The lower limit corresponds to the 
case of separation where the boundary layer equations are 
no longer valid. For the upper limit, the integral techni- 
que predicts velocity profiles which are unreasonable from 
a physical point of view. Therefore the solution is not 
valid in regions of highly adverse and highly favorable 
pressure gradients. The boundary layer solution is valid 
for laminar boundary layers in the presence of a laminar 
free stream. Although a laminar boundary layer can arise 
from the impingement of a turbulent jet, the free stream 

is turbulent which increases the heat transfer on the 
impingement surface. 

The finite element method proved very versatile in 
handling all the different types of jet flows considered. 
Although an iterative procedure was required in locating 
the free surface, the finite element method was efficient 
for Seiad these potential flow problems. Because of the 
complex nature of the jet impingement melting problem the 
integral solution provides the easiest method of calculating 
the boundary layer. Although the accuracy is affected by 
the approximate nature of the integral technique, the 
averaging effect of integrating the boundary layer equations 


produces reasonable results. 
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The conclusions from this work will be summarized 
in point form. 

1. The finite element method is generally known to 
be an effective means of solving field problems with complex 
boundaries. In this thesis, it was shown that the finite 
element technique is also a convenient method to use when 
the boundaries are not fixed a priori but are determined by 
the flow field itself. This was the situation for the free 
surface and the melting surface. Care must be taken, how- 
ever, to ensure that there are a sufficient number of nodes 
to accurately represent the flow field. In general more 
nodes are required when there are large changes in the velo- 
city or when the shape of the boundary changes rapidly. 

2. For a stagnation point flow the heat transfer 
rate is directly related to the square root of the main- 
stream velocity gradient along the impingement surface. 
This correlation provides a means of qualitatively measur- 
ing the heat transfer rates on the impingement surface. 
For example, in the case of a jet impinging into a cavity 
it was observed that the velocity gradient at the stagnation 
point was significantly reduced when the radius of curvature 
of the cavity tip was of the order of the radius of the jet. 
A corresponding reduction in the rate of heat transfer 
occurs. Also in the case of a jet impinging on a flat sur- 
face for nozzle-plate spacings of one jet diameter or less 
significant increases in the velocity gradients occurred 


near the lip of the nozzle. Again a region of high heat 
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transfer exists near the lip of the nozzle. 

3. There are three effects that make the heat trans- 
fer rate at the stagnation point of a water jet melting a 
cavity in ice different from that which would normally be 
associated with stagnation point flow. First there is the 
effect of injection of cold water into the boundary layer 
that is. caused by the melting. It was. found that this effect 
produces a decrease in heat transfer rate of up to sixty 
percent for very hot water in the jet. The second effect 
is produced by the highly temperature dependent properties 
of water near the freezing point. Although this effect was 
not analysed in detail it appears that the choice of the 
reference temperature at which fluid properties are calcu- 
lated will result in a fifty percent uncertainty in the 
dimensionless heat transfer rates. The final effect is 
produced by the fact that the shape of the cavity melted 
in the ice affects the flow field at the stagnation point 
and thus the heat transfer rate. On the basis of the calcu- 
lated heat transfer to a parabolic shaped cavity a signifi- 
cant effect of the cavity shape may have been expected. 
However, the cavity formed by melting was found to have a 
very flat bottom and the resulting change in heat transfer 
rate was only about ten percent at most. 

4. When cutting frozen sand or gravel the effects of 
shear stress are probably important. The present work 
assumed direct proportionality between cutting rate and 


shear stress, although the constant of proportionality was 
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unknown. 

This work is a start to predicting the melting rates 
of water jets impinging on ice surfaces. This knowledge 
could have practical application wherever ice formation was 
a problem and its removal necessary for the continued opera- 
tion of an industrial machine or component. Further to 
this work, the effects of turbulence on the melting process 
should be investigated. Also experimental work should be 
carried out to determine: 

a. which reference temperature should be used for 

calculating the physical properties; 

b. the value of the constant of proportionality c, 

when concerned with shear stress as a removal 


mechanism. 
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APPENDIX 1 


FINITE ELEMENT MATRIX FORMULATION 


In the finite element method described in Chapter III 
the solution of the potential flow field depended on the 
formulation of element matrices. By the use of area co- 
ordinates, these matrices may be found once and for all. 

In this appendix the element matrices for a general triangle 
as shown below for both two-dimensional and axisymmetric 
flows is briefly discussed. The notation of Chapter III 


ismused. This discussion follows that of Chan Plt]. 
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A.1 Matrices for Two-Dimensional Flow 


From Chapter III the element matrices are 
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The array ip is found by replacing the b's with a's in the 
expression for T,. 

Now substitution of the appropriate quantities into 
the matrix equations A.1 and A.2 enables their evaluation. 
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By taking the partial derivatives of the above integral 


with respect to dos $3 and oe results in 


] 
oo 
| on dy 
on 6 


S- 
fe 
Qa 
wn 
i 


om 
9¢. 


Q 
©6-|@ 
w 
eee | a ae | 
i 
I 
©- 
fob) Kod) 
= 
—J 
Q 
7) 
iT] 
Q 
=| 
[o>) 
—! 


Similar results can be obtained by considering the other 


two integrals. In this way the load matrix is obtained as 


= do 2 06,3 
SL, as ( (=) Lo 7 (a) L,)/6 


sty = ((2¢)32, + (38)'2,)/6 
SL, = ((22) 12, + (22)79,)/6 oad 
Sly = 2(22)%2,/3 


Par 


Isipasmt avods edd, tor as ae 
nt adfueey gb DAB pd og 


sotto of% paivsb tenes wit ioe om : 
26 bantetde eainaieviadniiine tr 


rit plat, 4 


weit) ) 


162 


Sher 2(22)'2/3 


”n 
rc 
W 


dpy2 
at) O5/3 


A physical interpretation of the load matrix may be obtained 


by considering the flux of velocity across a boundary surface. 


A.2 Matrices for Axisymmetric Flow 
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r= rye) + robo - r363 A.14 


The procedure is exactly the same as for the two-dimensional 
case and the element matrix and load matrix for axisymmetric 
flows are listed below. Here Pai is used to represent the 
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APPENDIX 2 


COMPUTER PROGRAM AND CUBIC SPLINE 
INTERPOLATION 


In this appendix the full computer program used to 
calculate the potential flow of an impinging jet and the 
resulting boundary layer is presented. A flow chart and 
a description of some of the important subroutines is given 
to aide in the understanding of this program. Also given 
is a description of the cubic spline interpolation technique 
used to calculate velocities and gradients of velocities on 
the impingement surface. The flow chart for the computer 
program is given in Figure A.1 and a description of the 
flow chart follows. 

The first step is to read the input data required for 
the finite element and boundary layer analysis. As dis- 
cussed in Chapter III only the nodal point data of the first 
six elements is required. Subroutine NGEN calculates the 
element node matrix for the remaining elements. Also formed 
are the matrices which consist of the nodes on various 
parts of the boundary. The X and Y or R and Z coordinates 
of the nodes on the boundary are read in and the coordinates 
of all other nodes are calculated in subroutine COORD. 

Every time a boundary node is moved, as in the adjustment 


of the free surface subroutine COORD must be called to recal- 
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Read input 
data 
NGEN 


Calculate coordinates 
of all nodal points 
COORD 


Generate element matrices 
MATGEN+LOAD 
Boundary conditions BC34 


Solve system of linear 
banded symmetric equations 
BSSOLV 


Calculate velocities 
on the free surface 
VELG 


Is free surface 
boundary condition 
satisfied? 


Adjust position of 
free surface 
ADJUST 


Calculate velocity on 
impingement surface 
UX 


Calculate the 
boundary layer 
MELT 


Is iteration 
concluded? 


MOVE 


Figure A.l Computer program flowchart 


vet 


to mebtieng feupba 
soattes e9%t 
TOULGA 


Saas 


. ore erry vr 


p bus p onan u 


Lt 


ae | 
a 2 - 
’ . 


at 4 uae 
; Weng, ui? bedi 
via 
Te re ) 


a) did usevll: a 
_ tasdowo St semowt setegmed, SoA eaigtt yg 
: a” vA f ‘y : 


168 


culate the coordinates of all the other nodes. 

Next the element matrices are calculated in subroutines 
MATGEN and LOAD while special boundary conditions are gene- 
rated in BC3. The resulting set of linear, symmetric banded 
equations are solved by subroutine BSSOLV. Depending on the 
form chosen the output from BSSOLV is either the nodal point 
velocity potential or the nodal point stream function value. 
The velocities at each of the nodes are now calculated in 
subroutine VELG. Subroutine ADJUST checks the velocities 
on the free surface with the velocity at the lip of the 
nozzle and adjusts the coordinates of the free surface node 
in the manner described in Chapter III. If the aves surface 
boundary condition is satisfied the program carries on, if 
not then with the adjusted free surface the program recal- 
culates the potential flow. This iteration procedure for 
finding the free surface is continued until the free surface 
boundary condition is satisfied. 

Subroutine UX then obtains the tangential and normal 
velocity components at points on the impingement surface 
from the X and Y or R and Z components of velocity calcu- 
lated by sipeeacene VELG. This in preparation for the bound- 
ary layer calculation contained in subroutine MELT. In this 
subroutine the position of the melted surface is calculated 
and subroutine MOVE transforms the melted surface as des- 
cribed in Chapter VII. Now one iteration of the program 
has been completed. The option for continuing the melting 


process is available, and the computer program proceeds to 


. geben ysdto 
zentivetdue nat borstuolss sis. 
-snep o76 enotbttbnos X14, bavod 
bebnsd Iiitemaye .1S63Ntt Yo aids 17 
ait oe garbneqed Vi0228 Raab ving eo 
tn toq febon oti vedtts et vuoee8 most 
.ouley nofsonet meats Intog {bon ode 7 
ant botefuolso woe 946 zebon ots to dace 26 
eotitoolev ead agoedo “TBULGA sotaverdua . D4aV an 
ait to ght oat f6 vttsotew. odd. ddiw gostave wit 
Shon sostiue S877 sit to ‘setbotarege edt cesarean 
sonttive sont Ad tI tt vedigsrtd at: poatroasd 19 
+h to ee hrred ae719019 of? bettatiee et sorsione9 x "es 
-fso94 aerporg ai? sosttue ‘99%? bos eutbs ods Ad tw 
sot sr ebssarg notiensat erat ere Fettassog Preae 
gostu2 seth ont, | tomu bavarinos ci ‘were ae1t sae) 
ihettettes « eF notsthaoa ¥ 
fsarrvon bas [ettnepaes dikes ois XU entiuor 
aseTAve Sevgamsigeet yt aan tie azatoa “Ye etnenogmos arcolon' | 
-ualeo ¢tfoolsy, Yo <3ne | ie RCH eT RE seg 
-bauod odt tot Hettersqa4g AF .019¥ satinordue Yd base } 
2tdd ol .TIaM enttuoydue wi bantstno. vottetvates 13ysl ts | 
bsiefuaten 2f sostrse badttsm os Yo notireoq od? ‘ontauordue | rh 


-# 


ze. 


-29b) és) aos7av2 bottom sad sovotencys 3VOM enttuorduz bas 
msipeyq SAF To hott siatit ga woh 114. as el nt hatte 
oni d fom afd onfudtrs nos 10% potzqo sat beat gmoo "naod 26a 


a3 abasn0%d maypord Metuaiios on9 bos .sidetteve et 2292070 


169 


calculating the matrices for the potential flow on the 
adjusted melted surface and the process is repeated. 

In order to understand the boundary layer calculation 
a brief description of subroutine MELT is given. The input 
to this subroutine are the coordinates of the points on the 
boundary and their tangential and normal velocities. First 
the coefficients required in the cubic spline interpolation 
of the velocities are calculated in subroutine BCS. The 
cubic spline interpolation routine is described in some 
detail by Ahlberg, Nielson and Walsh [55]. The parameters 
at the stagnation point are then evaluated. 

Between each pair of nodal points the boundary layer 
parameters are calculated at ten equally spaced intervals. 
The velocity and gradient of velocity at these interior 
points are calculated by subroutine SPFUN using the cubic 
spline coefficients from subroutine BCS. Subroutine ROOT 
calculates the value of the new pressure gradient parameter 
A at each successive point on the boundary layer. The posi- 
tion of the melted surface is calculated at the nodes at 
the boundary as the solution procedure marches along the 
impingement surface. 


The computer program follows. 
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JET IMPINGEMENT HEAT TRANSFER PROGRAM 
POTENTIAL FLOW ANALYSIS COMPUTED BY 
THE FINITE ELEMENT METHOD 

BOUNDARY LAYER ANALYSIS COMPUTED BY 
THE KARMAN POHLHAUSEN INTEGRAL METHOD 


ie 3 FE A I a a aie ae i fe 2 SE ae 2K a I a ic 9 3 2k le 2 9 i a i 2 2 a aK a 2 a 2k i 5 3 I IK 2 CK CE OK 2K 2K 


MAGA ath aaa a 


REAL S (3,3),SA (6,6) ,SL(6) -7 (6,6) ,TH (6,6) ,SUT(3) ,SUN(3), 
*VELX(6) ,VELY (6) , VEL(6) , THETA (6) ,SIDE(3) ,CHI(6,3) , DX (3), 
*X (343) ,¥ (343) ,GS (343,17) ,GL (343) , VI (27) , VN(27) ,SX (27) 
*,DY(3) ,SZ(27) ,SR(27) ,H(81) ,GG (27) ,DG (27) ,ZCOR (27) 

INTEGER LEL(144,6) , NAB(23),NBC(27) ,NDE(60) ,NEF (49), 

*N4 (25,4) ,NCD (7) ,NFA(7) ,NBE(3) , NBN (7) 

COMMON /ELE/ L,DX,DY,Y1I,Yd,YK,X1,XJ,XK, AREA 

COMMON /SPEED/ CHI, VELX, VELY,VEL,THETA,T,TH,NB 

COMMON /BOUND/ KS,KE,XST 


READ INPUT DATA 

NE=THE NUMBER OF ELEMENTS 

NN=THE NUMBER OF NODES 

NA=G TWO DIMENSIONAL FLOW 

NA=1 AXISYMMETRIC FLOW 

NB=G STREAMFUNCTION FORMULATION 

NB=1 POTENTIAL FUNCTION FORMULATION 
NI=THE NUMBER OF ITERATIONS FOR FREE SURFACE 
KS=THE STAGNATION NODE NUMBER 

KE=THE LIP OF THE NOZZLE NODE NUMBER 
KBW=THE EXPECTED BANDWIDTH 

NT=THE NUMBER OF MELTING ITERATIONS 
NDIV=THE NUMBER OF POINTS BETWEEN NODES 
NSS=0 NO SHEAR STRESS 

NSS=1 SHEAR STRESS ACTING 


Goeaeadgsaanannavraanacens 


READ (5,900) NBD,NA,NB,NI,KS,KE 

READ(5,900) KBW,NT,NDIV,NSS 
900 FORMAT (6110) 

IER=0 

NE=6* NBD 

NN=7+14*NBD 

NNG=14+NBD 

NABN=1+(KS-1) /7 

NBCN= (NN-6-KS) /7+1 

NB3=3*NBCN 

NCDN=7 
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NDEN= (NN-KE) /7 

NFS=NDEN/2 

NEFN= (KE-7) /7+1 

NFAN=7 

NBEE=3 

NBNN=7 

READ(7,902) AMDA,EPS,XST,VO 
READ(7,902) PR,STE,21,GAMMA 
READ(7,902) DT 

FORMAT (4F20. 6) 

READ(7,901) ((CHI (I,J) ,J=1,3) ,I=1,6) 
FORMAT (3F20.6) © 

READ(7,900) ({(LEL (I,J) ,J=1,6) ,I=1,6) 
CALL NGEN({NE,LEL,NNG, N4, NN,X,Y,NABN,NAB,NBCN,NBC,NCDN, 
*NDE,NDEN, NDE, NEFN, NEF,NFAN, NFA, NBE, NBN) 
DO 10 I=1,NABN,2 

NO=NAB (I) 

READ (7,901) X (NO) 

CONTINUE 

DO 11 I=1,NBCN,2 

NO=NBC (I) 

READ (7,901) X (NO) ,¥ (NO) 

CONTINUE 

DO 12 I=2,NDEN,2 

NO=NDE (1) 

READ(5,901) X (NO) ,Y¥ (NO) 

CONTINUE 

DO 14 I=1,NEFN,2 

NO=NEF (I) 

READ(5,901) X(NO) 

CONTINUE 

CALL COORD(N4,NN4,NN,NE,LEL, Y, X) 
CALL BAND (NBW, NN, NE, LEL) 

IF (NBW.GTI.KBW) GO TO 300 


START MELTING ITERATIONS 


pO 500 IJK=1,NT 


START ITERATIONS TO FIND FREE SURFACE 


DO 1 LL=1,NI 
LNI=LL 


CALCULATE ELEMENT MATRICES 
INSERT BOUNDARY CONDITIONS 


DO 17 M=1,NN 
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DO 17 N=1,NBW 

GS (M,N) =0.0 

DO 15 M=1,NN 

GL (M) =0.0 

LMM=0 

DO 20 L=1,NE 

CALL NODE(NN,NE,LEL,Y,X, IER) 
IF(IER.EQ.1) GO TO 301 

CALL MATGEN({GS,NN, NE, LEL, NA, NBW) 
DO 20 I=1,NBEE 

IF (L.EQ.NBE(I)) GO TO 22 

GO TO 20 

CALL LOAD(GL,NN,NE,LEL, NA, NBN, NBNN, VO) 
CONTINUE 

CALL BC3(GS,GL,NCD, NCDN, NN, NBW) 


SOLVE LINEAR BANDED SYMMETRIC SYSTEM 


CALL BSSOLV(NN,NBW,GS, GL) 


CALCULATE NODAL POINT VELOCITIES 


DO 29 I=1,NN 

DO 29 J=1,4 

GS (I,J) =0.0 

IF(J.EQ.4) GS(I,J) =1.0 
CONTINUE 

DO 30 L=1,NE 

CALL NODE(NN,NE,LEL,Y,X, LER) 
IF(IER.EQ.1) GO TO 301 

CALL VELG(NN,NBW, NE, LEL, GL, GS) 
CONTINUE 

KE6=KE-6 

VIA=GS (KE, 3) 

UF=VIA/GS (KE6,3) 

DO 60 M=1,NN 

DO 60 124,3 

GS (M,I)=GS(M,1I)/VIA 

VIA=1.0 


ADJUST FREE SURFACE 


CALL ADJUST (NDE, NDEN, VIA, Y,X,GS,NN,NBW,LMM,AMDA, EPS) 


IF(LMM.EQ.NFS) GO TO 100 
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CALL COORD(N4,NNU, NN, NE, LEL, Y,X) 

CONTINUE 

GO TO 302 

WRITE (6,800) LNI 

FORMAT ('O',* FREE SURFACE LOCATED # OF ITERATIONS=',14) 
WRITE (6,904) UF 

FORMAT('O',* FREE SURFACE VELOCITY',F10.6) 

WRITE (6,801) 

FORMAT (*O0",*NODE',11X," VELOCITY! ,20K,'X",22Xx,'Y') 
DO 70 J=2,NDEN,2 

LJ=NDE(J) 

WRITE (6,802) LJ,GS(LJ,3) ,X(LJ) ,¥ (1d) 

FORMAT (* *,13,5X,F15.6,9X,F1506,9X,F15.6) 


PREPARE VELOCITIES FOR BOUNDARY LAYER 
CALL UX (NBC, NBCN,X,Y,SX,SZ,S8,GS,NN,NBW,VT,VN,VIA) 


CALCULATE BOUNDARY LAYER 


CALL MELT (NBCN,VT,SX,S2,SR,PR,STE,R1,GAMNA,NA,NDIV, 


*ZCOR, IER,UF, DT,NSS) 


IF (IER.GT.0) GO TO 999 

WRITE (6, 198) 

FORMAT ('O',* POSITION OF MELTED SURFACE') 

WRITE (6,199) (NBC (I) ,ZCOR (I) ,¥ (NBC(I)) ,1=1,NBCN) 
FORMAT (*0',3 (15, 2F15.6)) 

CALL MOVE (ZCOR,NBC, NBCN,X,Y, NN, NDE, NDEN) 

CALL COORD(N4,NN4,NN,NE,LEL, Y,X) 

WRITE (6,200) 

FORMAT (*0',* RELOCATED COORDINATES ON THE SURFACE') 
WRITE (6,201) (NBC (I),X (NBC(I)),¥ (NBC(1)) , 1=1,NBCN) 
FORMAT (*0',3 (15, 2F15.6)) 

CONTINUE 


GO TO 999 


WRITE (6,400) NBW 

FORMAT('1",* BANDWIDTH TOO LARGE',I4) 

GO TO 999 

WRITE (6,401) IER 

FORMAT ('1",*NEGATIVE AREA IER=', 13) 

GO TO 999 

WRITE (6,402) LNI 

FORMAT('1'," FREE SURFACE NOT FOUND WITHIN',13,2X, 


**TTERATICNS') 


CONTINUE 
po 39 I=1,NBCN,2 
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NO=NBC (I) 

WRITE (10,901) X(NO),¥ (NO) 
DO 40 I=2,NDEN,2 

NI=NDE (I) 

WRITE (10,901) X(NI) ,Y(NI) 
STOP 

END 
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(C4 2 3g 9 a aie 2 a a 2 a a 2c i 2 ic aK 2 9 ig aie aK aie ie a 2 a 2 i ak aK a Tig DC IC ie Ak RAI IE 3c ik 38 2c ik a 2k 2c iC OK OK 3 
Cc 
3 SUBROUTINES FOR FINITE ELEMENT PROGRAM 
Cc 
(CE HR A I a 2K 2 aK 2 a 2 i 2 a eI i 2 iC I a aK aC 9K 2 9 a 2 Ak 2 He IK TE IK 3K 3k KK 2k 2 2k ak a OK Ag 2 2c OK a 
SUBROUTINE CCORD (NU, NNU, NN, NE, LEL,Y,X) 
REAL ¥ (NN) ,X (NN) 
INTEGER N4(NN4,4) , LEL (NE, 6) 
DO 300 N=1,NN4 
X (N4 (N, 3)) =(X (N4 (N, 1) ) +3. O*X (N4(N,4))) #0. 25 
¥ (N4 (N,3)) = (CY (94 (MN, 1) ) +3. 0*Y (N4 (N,4))) #9. 25 
¥ (NY (N, 2)) =(3« O*Y¥ (NG (N,1)) +¥ (NG (N,4))) *0. 25 
X (NG (N, 2) ) =(3-0*K (N4 (N,1)) #X (NY (N, 4) )) #0. 25 
300 CONTINUE 
DO 301 L=1,NE 
DO 301 I=1,3 
J=iI+1 
if (J.,EQ. 4) J=1 
K=I+3 
LI=LEL (L,1) 
LJ=LEL (L,J) 
LK=LEL (1, K) 
Y (LK) =0.5%* (¥ (LI) +¥ (LJ) ) 
X (LK) =0. 5* (X (LI) +X (LJ) ) 
301 CONTINUE 
RETURN 
END 


SUBROUTINE BAND(NBW, NN, NE,LEL) 

INTEGER LEL(NE,6) 

NBW=1 

DO 1-L=1,8E 

poi) T=1;6 

DO 1 J=1,5 

LI=LEL (L,1) 

LJ=LEL(L,J) © 

IF (LI.LT.1J) GO TO 1 

NCOL=LI-LJ+1 

IF (NBW-NCOL.1T.0) NBW=NCOL 
1 CONTINUE 

RETURN 

END 


SUBROUTINE NODE(NN, NE,LEL,Y,X,1ER) 

REAL Y(NN) ,X (NN) ,DY(3) ,DX (3) 

INTEGER LEL(NE,6) 

COMMON /ELE/ L,DX,DY,YI,YJ,YK,XI,XJ,XK,AREA 
YI=Y (LEL (L,1)) 
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¥I=¥ (LEL@L,2)) 

YK=¥ (LEL (L,3)) 

XI=X (LEL (£1, 1)) 

XJ=¥(LEL (1 2)) 

XK=X (LEL (L,3)) 

DX (1) =XK-XJ 

DX (2) =XI-XK 

DX (3) =XJ-XI 

DY (1) =YJ-YK 

DY (2) =YK-YI 

DY (3) =YI-YJ 

AREA= (DX (3) ¥ DY (2) -DY (3) * DX (2)) /2.0 
IF(AREA.LE.0.0) GO TO 7 

RETURN 

WRITE(6,100) L,AREA 

AREA=-AREA 

IER=1 

FORMAT ('O*", "ELEMENT',110,' NEG AREA',F 16.3) 
RETURN 

END 


SUBROUTINE MATGEN(GS, NN, NE, LEL, NA, NBW) 
REAL S (3,3) ,SA (6,6) ,DX(3) ,DY (3) ,GS (NN, NBW) 

INTEGER LEL(NE,6) 

COMMON /ELE/ L,DX,DY,Y1,YJ,YK,X1,XJ3,XK,AREA 

RI=YI 

RJ=YJ 

RK=YK 

IF(NA.5BQ.0) GO TO 16 

DO 13 “=1,3 

DO 13 J=1,3 

5 (I,J) = (DX (I) *DX (J) #DY (I) *DY (J) ) / (60. 6* AREA) 

SA (1, 1) =3*S (1,1) * (3*RI+RJ+RK) 

SA (1, 2) =-S (1,2) * (2* RI +2* RJ+ BK) 

SA (1,3) =-S (1, 3) * (2* RI + RI +2*RK) 

SA (1,4) =S(1, 1) *(3*RI-2*RJ-RK) +S (1,2) * (14*RI+3*RI+3*RK) 
SA (1,5) =S (1, 2) ¥ (3¥ RI-RJ-2*BK) +5 (1, 3) * (3*¥RI- 2*RI-RBK) 

SA (1,6) =S (1, 1) * (3*RI-Rd-2*RK) +S (1,3) * (14*RI+ 3*RJ+3*BK) 
SA (2, 2) =3*S (2,2) * (RI+3*RJ+RK) 

SA (2,3) =-S (2,3) * (RI +2* RJ +2* BR) 

SA (2,4) =S(1, 2) ®(3*RI+14*RI+3¥*RK) +S (2,2) * (-2*RI+3*RJ-RK) 
SA (2,5) =S (2, 2) ¥(—RI+3*RJ -2*RK) +5 (2, 3) * (3*RI414*RJ+ 3* RK) 
SA (2,6) =S (1, 2) ¥ (-RI+3*RJ-2*RK) +S (2, 3) * (-2*RI + 3*RJ-RK) 
SA (3, 3) =3*S (3,3) * (RI+RJ+ 3*RK) 

SA (3,4) =S (1, 3) ¥ (-RI-2* RJ +3*RK) +5 (2, 3) * (- 2#RI-RI+3*RK) 
SA (3,5) =S (2, 3) * (3¥RI+ 3*RI+14*RK) +5 (3,3) * (-RI-2¥*RJ+ 3*BK) 
SA (3,6) =S (1, 3) ¥(3*RI+3*RI+14*RK) +S (3, 3) * (-2*RI-RI+3*RK) 
SA (4,4) =8* (S (1,1) * (RI+3*RJ+ RK) -SA(1, 2) +S (2, 2) * 
* (3*RI+RJ +RK) ) 
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SA (4,5) =8*S (1,3) * (RI+3*RI+RK) —4* (S (1, 2) ¥RI+ 
*S (2,2) *RJ+S (2,3) * BK) 

SA (4,6) =8*S (2, 3) * (3¥RI+RJ+RK) -4* (S (1,1) *RI+ 
*S (1,2) *RJ+S (1,3) *RK) 

SA (5,5) =8*(S (2,2) * (RI+ RJ +3*RK) -SA (2,3) +S (3,3) * 
* (RI+3*RJ+RK) ) 

SA (5,6) =8*S (1,2) * (RI+RJ+3¥*RK) —4* (S (1,3) *RI+ 
*S (2,3) *RJ+S (3,3) *BK) 

SA (6,6) =8* (S (1,1) * (RI+RJ4+3*RK) -SA (1, 3) +S (3, 3) * 
* (3*RI+RJ+RK) ) 

Go T0 17 

CONTINUE 

DO 20 -1=1,3 

pO. 20 “3773 

S (1, J) = (DX (1) * DX (J) +DY (I) *DY (J) ) / (12. O* AREA) 
SA (1, 1) =3.0*S (1,1) 

$A (4; 2)]=5 (1,2) 

SA (1,3)=-S (1,3) 

SA (1,4) =4.0*S (1,2) 

SA (1,5) =0.0 

SA (1,6) =4.0*5 (1, 3) 

SA (2,2) =3.0*S (2,2) 

SA (2,3) =-S(2,3) 

SA (2,4) =4.0*5 (1,2) 

SA (2, 5) =4.0*S (2,3) 

SA (2,6) =0.0 

SA (3, 3) =3.0* 5 (3,3) 

SA (3,4) =0.0 

SA (3,5) =4.0*5 (2, 3) 

SA (3,6) =4.0*S (1,3) 

SA (4, 4) =8.0* (S (3,3) -S (1,2)) 

SA (4,5) =8.0*S (1,3) 

SA (4, 6) =8.0*S (2, 3) 

SA (5,5) =8.0* (S (1,1) -S (2, 3)) 

SA (5,6) =8.0*S (1,2) 

SA (6,6) =8.0* (S (2,2) -S (1, 3)) 

CONTINUE 

pO 18 mM=1,6 

pO 14 N=1,6 

SA (N,M)=SA(M,N) 

po 11 I=1,6 

pO 11. J=1,6 

LI=ERU GL (1) 

LJ=LEL (L, J) 

IF(LJ.LT.L1I) GO TO 11 

K=LJ-LI+1 

GS (LI, K) =GS(L1I,K) +SA (I,J) 

CONTINUE 

RETURN 

END 
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SUBROUTINE LCAD(GL,NN,NE,LEL,NA,NBN, NBNN,VO) 
REAL GL(NN) , DX (3) , DY (3) , SIDE (3) ,SL (6) 
INTEGER LEL(NE,6) , NBN (NBNN) 

COMMON /ELE/ L,DX,DY,Y1,YJ,Y¥K,X1,XJ,XK,AREA 
DO 38 I=1,3 

S IDE (I) =SORT (DX (I) **2+DY (I) **2) 
SIDE (1) =0.0 

SIDE (3)=0.0 

IF (NA.EQ.0) GO TO 40 

SL (1) =YI*VO* (SIDE (2) +SIDE (3) ) /6.0 
SL (2) =¥d *VO* (SIDE (3) +SIDE (1)) /6.0 
SL (3) =YK*VO* (SIDE(1) +SIDE (2) ) /6.0 
SL (4) = (YI+YJ) *VO*SIDE (3) 73.0 

SL (5) = (Yd+¥K) *VO*SIDE (1) /3.0 

SL (6) = (YI+YK) *VO*SIDE (2) /3.0 

GO TO 41 

SL (1) =VO* (SIDE (2) +SIDE(3)) /6.0 

SL (2) =VO* (SIDE (3) +SIDE(1)) /6.0 

SL (3) =VO* (SIDE (1) +SIDE(2)) /6.0 

SL (4) =2. O¥VO*SIDE (3) /3.9 

SL (5) =2. O¥VO*SIDE (1) /3.0 

SL (6) =2.0*V0 ¥SIDE (2) /3.9 

CONTINUE 

DO 42 M=1,6 

IF (L.LT.20) SL(M) =-SL(M) 

CONTINUE 

DO 43 I=1,6 

DO 43 N=1,NBNN 

LI=LEL(L,1) 

NBB=NBN (N) 

IF (LI.EC.NBB) GL(NBB) =GL (NBB) +SL (1) 
CONTINUE 

RETURN 

END 


SUBROUTINE BSSOLV (NN, NBW, GS, GL) 
REAL ST(18),GS(NN,NBW) ,GL (NN) 
N=9 
N=N+1 
REDUCE PIVOT EQUATION 
GL (N) =GL (N) /GS (N,1) 
IF (N-NN) 150,309,150 
DO 200 K=2,NBW 
ST (K) =GS (N,K) 
GS (N,K) =GS(N,K) /GS (N,1) 
REDUCE REMAINING EQUATIONS WITHIN SPAN 
DO 260 L=2,NBW 
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I=N+L-1 
IF(NN-I) 260,240,240 
J=9 
DO 250 K=L,NBW 
J=I3+1 
GS (I,J) =GS{1,J)-ST({L) *GS (N, K) 
GL (I) =GL (1) -ST (L) *GL(N) 
CONTINUE 
GO TO 100 
BACK SUBSTITUTION 
N=N-1 
IF(N) 350,500,350 
DO 400 K=2,NBW 
L=N+K~-1 
IF (NN-L) 400,370,370 
GL (N) =GL (N) -GS (N, K) *GL(L) 
CONTINUE 
GO TO 300 
CONTINUE 
WRITE (10,900) 
FORMAT ('-',* CALCULATED NODAL VALUES’) 
DO 600 N=1,NN 
WRITE (10,901) N,GL(N) 
FORMAT ('0',110,F20.6) 
CONTINUE 
RETURN 
END 


SUBROUTINE VELG(NN,NBW,NE,LEL,GL,GS) 

REAL VELX (6) ,VELY (6) , VEL (6) , THETA(6) , DX (3) ,DY¥ (3), 
*T (6,6) ,TH (6,6) ,GS(NN,NBW) ,GL(NN) ,CHI (6, 3) 
INTEGER LEL(NE,6) 

COMMON /SPEED/ CHI, VELX,VELY,VEL,THSTA,T,TH,NB 
COMMON /ELE/ L,DX,DY,YI,YJ,YK,X1,XJ,XK, AREA 
PI=ARCOS {-1. EO) 

pO 2 I=1,6 

DO 2 J=1,3 

T (I,J) = (4.0*CHI (I,3)-1.0) *DY (J) / (2. 0O* AREA) 

TH (1,J)= (4.0 *CHI (I,J) -1. 9) ¥DX (J) / (2. O* AREA) 
CONTINUE 

DO 3 I=1,6 

DO 3 J=1,3 

JJ=J+3 

K=J+1 

IF (Ke-EQ.4) K=1 

T(I,J3J) =2.0* (DY (J) *CHI (I,K) +DY (K) *CHi (I,J) ) /AREA 
TH (1, Jd) =2.0 * (DX (J) *CHI (I,K) +DX (K) *CHT(1,J))/AREA 
CONTINUE 

po 4 I=1,6 
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SUMX=0.0 

SUMY=0. 0 

DO 5 J=1,6 

LISLEL (1,0) 
SUMX=SUMX+4GL (LJ) *T (I,J) 
SUMY=SUMY+GL (LJ) *TH (I,J) 
CONTINUE 

VELX (I) =SUMX 

VELY (1) =SUMY 

IF (NB.EQ.0) VELX (I) =SUMY 
IF (NB.EQ.0) VELY (I) =-SUMX 
VEL (I) =SQRT ( VELX (I) **2+VELY (1) **2) 
THETA (I) =180.0*PI 

LI=LEL (L,I) 

IF (GS (11,3). 8E.0-0) GO TO 11 
GS (LI, 1) =VELX (1) 

GS (LI, 2) =VEL Y (I) 

GS (LI,3)=VEL (1) 

GO TO 4 

GS (LI, 4) =GS (11,4) +1.0 
VX=VELX (I) 

VY=VELY (I) 

VF=VEL (1) 


GS (LI, 1) = (VX+GS (LI, 1) * (GS (LI, 4) -1.0)) /GS (LI, 4) 
GS (LI, 2) = (VY+GS (LI, 2) * (GS (LI,4) -1.0)) /GS (LI,4) 
GS (LI, 3) = (VF+GS(LI, 3) * (GS (LI,4)-1.0) )/GS (LI,4) 


CONTINUE 
RETURN 
END 
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SUBROUTINE NGEN(NE,LEL,NN4,N4Y,NN,X,Y,NABN,NAB,NBCN,NBC, 


*NCDN, NCD, NDEN, NDE, NEFN, NEF, NFAN, NFA, NBE, NBN) 


REAL X (NN) ,Y (NN) 


INTEGER LEL(NE,6),N4(NN4,4) ,NBE(3),NBN(7) ,NAB(NABN), 
*NBC (NBCN) ,NCD (NCDN) ,NDE(NDEN) , NEF (NEFN) , NFA(NFAN) 


COMMON /BOUND/ KS,KE,XST 
NBEN=3 

NBNN=7 

DO 2 L=7,5E 

J=L-6 

po 2 11,6 

LJ=LEL (J,1) 
LEL(L,1)=Ld+14 

pO 3 I=1,NN4 

J= (1-1) * 14-1 


po 3 K=1,4 
J=I5+2 
N4 (I,K) =d 


K=NE-6 
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po 4 I=1,3 
II=1*2 

NBE (I) =I1 
K=NN-7 

po 5 i=1,7 
NFA (I) =1 

NBN (1) =1 
K=K+1 
NCD(1) =K 

DO 6 I=1,NABN 
K= (1-1) *7+1 

¥Y (K) =0.0 
NAB(I) =K 

DO 7 I=1,NBCN 
K= (I-1) *7#KS 
NBC (I) =K 

pO 8 I=1,NDEN 
K=KE+7*1 

NDE (I) =K 

DO 9 I=1,NEFN 
K=7+7* (I-1) 

¥ (K) =0.5 

NEF (I) =K 
WRITE (10,800) 
FORMAT (6115) 
WRITE (10,800) 
WRITE (10,809) 
WRITE (10,800) 
WRITE (10,800) 
WRITE (10,800) 
WRITE (10,800) 
WRITE (10,800) 
WRITE (10,800) 
WRITE (10,800) 
RETURN 

END 


SUBROUTINE ADJUST(NDE,NDEN,VIA,Y,X,GS,NN, NBW,LMM,AMDA 
* , EPS) 


REAL Y(NN) -X ( 
INTEGER NDE(N 
DO 1 J=2,NDEN 
II=J-1 

K=J+1 

IF (J.EQ» NDEN) 
LI=NDE (11) 
LJ=NDE (J) 
LK=NDE (K) 


WRITE (8,101) LJ,GS(LJ,3) ,VIA 


((LEL(1I,J) ,J=1,6) , 1=1, NE) 


((N4 (I,J) , d=1,4) ,1=1,NNA4) 
(NBE (I) ,I=1,NBEN) 
(NBN(I) ,I=1,NBNN) 
(NAB (I) ,I=1,NABN) 
(NBC (I) , I=1,NBCN) 
(NCD (I) ,I=1,NCDN) 
(NDE (1) ,1=1,NDEN) 
(NEF (I) ,I=1,NEFN) 
(NFA (I) ,I=1,NFAN) 


NN) ,GS(NN,NBW) 
DEN) 
22 


K=J 
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101 FORMAT (*'O','"NODE',I5,'VELOCITY',F12.5,'ACTUAL VELOCITY? 
*,F12.5) 
¥II=Y (11) 
Y3d=Y (LJ) 
YK K=Y (LK) 
XII=X (LI) 
XJ J=X (LJ) 
XKK=X (LK) 
WRITE (8,102) LdJ,XJJ,YJJ 
102 FORMAT(110,2F20.8) 
IF (YII-YKK.EQ.0.0) GO TO 2 
SLOPE= (XII-XKK) / (YII-YKK) 
GO TO 3 
2 SLOPE=1.05810 
3 CONTINUE 
IF (SLOPE.EQ.0.0) SLOPE=1.0E-10 
SLN=-1.0/SLOPE 
IF (ABS (GS(LJ,3)-VIA).LT.EPS) GO TO 10 
DIST=AMDA* ({GS (LJ, 3) /VIA) **2-1.0) 
D=DIST/SQRT (1. 0+SLN**2) 
IF (SLN.GE.0.0) GO To 9 
¥ (LJ) =Y (LJ) +D 
X (LJ) =X (1J) +D*SLN 
WRITE (8,102) LJ,X(LJ) ,¥ (LJ) 
GO 70 1 
9 CONTINUE 
Y (LJ) =¥ (IJ) -D 
X (LJ) =X (LJ) - D*SLN 
WRITE (8,102) LJd,X (LJ) ,¥ (LJ) 
GO TO 1 
10 WRITE (8,100) LJ,GS (LJ, 3) 
100 FORMAT('-','NO ADJUSTMENT NODE='!,15,"VELOCITY=',F12,5) 


LMM=LMM+1 
1 CONTINUE 
RETURN 
END 
c 
c 
SUBROUTINE BC3(GS,GL, NCD, NCDN,NN, NBW) 
. POTENTIAL FLOW BOUNDARY CONDITIONS 
Cc UNIFORM FLOW AT THE OUTLET 
c SPECIFIED VELOCITY AT THE INPUT FROM SUBROUTINE LOAD 


REAL GS (NN, NBW),GL (NN) 
INTEGER NCD(NCDN) 
| pO 1 I=1,NCDN 
K=NCD (1) 
K1=K 
GL {K) =0.0 
GS (K,1)=1.0 
DO 1 J=2,NBW 
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K1=K1-1 

GS (K1,J) =0.0 
Ss: (k 5} =000 
RETURN 

END 


SUBROUTINE UX(NBC,NBCN,X,Y¥,SX,SZ,SR,GS,NN,NBW,VT,VN,VIA) 
REAL X(NN),Y (NN) ,GS(NN,NBW) ,VT(NBCN) ,VN (NBCN) , SUT (3), 


*SUN (3) ,SX (NBCN) ,5Z (NBCN) ,SR (NBCN) 


INTEGER NBC (NBCN) 

DC 1 N=1,NBCN 

VT (N) =0.90 

VN (N) =0.0 

DO 2 J=2,NBCN,2 

I=J-1 

K=J+1 

II=NBC (I) 

KK=NBC (K) 

XI=X (11) 

YI=¥ (IT) 

XK=X (KK) 

YK=Y¥ (KK) 

THE=ABS (ATAN ( (XK-X1I) / (YK-YT) )) 
DO 2 M=1,3 

JJ=J+M-2 

NJJ=NBC (JJ) 

SUT (M) =GS(NJJ, 2) *COS (THE) -GS (NJJ, 1) *SIN (THE) 
SUN (M) =- (GS(NJJ,2) *SIN (THE) +GS (NJJ,1) *COS (THE) ) 
IF(VT(JJ).EQ.0.0) GOTO 4 

VT (JJ) = (VE (JI) +SUT (M)) /2.0 

VN (JJ) = (VN (Jd) +SUN(M)) 72.0 

GO TO 2 

VT (JJ) =SUT (M) 

VN (JJ) =SUN(M) 

CONTINUE 

SX (1) =0.0 

SZ (1) =4.9 

SR (1) =0.0 

DO 5 K=2,NBCN 

KI=NBC (K) 

XK=X (KI) 

Y¥ K=Y (KI) 

J=K-1 

JI=NBC (J) 

YJ=Y (JI) 

X J=X (JI) 

SZ (K) =X (KI) 

SR (K) =Y (81) 

SX (K) =SX (J) + SQRT ( (XK- XJ) **2+ (YK- YJ) **2) 
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WRITE (6,100) 
100 FORMAT ('-',*NODE',21X,'S',11X,'TANGANTIAL VELOCITY',7X, 
*'NORMAL VELOCITY") 
DO 3 I=1,NBCN 
NI=NBC (I) 
WRITE (6,101) NI,SX(I),VT(I),VN(I) 
101 FORMAT(* *,14,10X,F15.6,2(8X,F15.6)) 
3 CONTINUE 
RETURN 
END 
c 
(8 Be 2 eK ik 2 2 Sic ic 3K 2c ik 2k 2g 2k ie aK ik 2k ae 2k is 2 3c kkk ak ek Ck ik fe aC ae ake oii a ale a ik aie a ie aK oe a IS ic ae Kk ak ok 
Cc 


c BCUNDARY LAYER SUBROUTINE MELT 

c AND ASSOCIATED SUBROUTINES 

S 

CAR AR I OI 32 I I a I a aa a I a 2 RI I ARK IO ak 2 aI a aK 2c ak ak 2 ake 2 2 2c 
Cc 

c 


SUBROUTINE MELT (NBCN,U,S,ZA,R,PR,ST,R1,GAMMA,NA,M, ZCOR, 
*IER,UF,DT,NSS) 

REAL S(NBCN) ,U (NBCN) ,H(123),R(NBCN) ,GG(41),DG (41), 
*ZA (NBCN) ,ZCOR (NBCN) 

po 10 I=1,NBCN 

GG (1) =5 (Z) 

U (I) =U (I) *UF 

10 ZCOR(I)=ZA (ZI) 

G=GAMMA 

IER=0 

PRP=SQRT (0. 254+5T/3.0)-0.5 

N=NBCN 

NOP=N 

P=PRP/PR 

PP=P/R1 

WRITE (6,200) ST,PR,P,R1 

200 FORMAT(*O','STEFAN NUMBER',F6.2,/,'PRANDTL NUMBER',FO.2, 

*/, "MELTING PARAMETER", F10.6,/,'BOUNDARY THICKNESS RATIO! 
*,F10.6) 

H(1)=0.0 

ZS=0.0 

DZ 1=0.0 

wS=0.0 

DW1=0.0 

CALL BCS(U,S,N,NOP, GG, DG, H, IER) 

JJ=2 

X=0.0 

CALL SPFUN(H,U,S,N,X,JJ,V,DV) 

CALL FUN(PR,ST,R1,G, FG, GGG, GK, GL, D1, D2, D3) 
Z0=GK/D¥ 

WO=D3*D3¥*R1* B1*Z0/ (D2*D2) 
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IF (NA.EQ.0) ZS=2Z0 

IF (NA.ECQ.0) WS=WO 

IF (20.LT.0.0.0R.WO.LT.0.0) GO TO 300 
THK2=SQRT (Z0) 

THK3=SQRT (WO) 

BL=THK2/D2 

BL F=THK3/D3 

VW=6.0*P/BLE 

VW1=6,0*P/ (R1*BL) 

DUUDY= (2.0+G/6.0) / (BL+BL* EP) 
DTDY=2. 0/ (BLE+BLE*PRP) 

CTSR=V*DUUDY 

WRITE (8,212) VW,VW1 

FORMAT ('O',*VW,VW1",2F20. 6) 
DWZ=VW*DT+CTSR*DT*NSS 

ZA (1) =ZA (1) +DWZ 

ZCOR (1)=ZA (1) 

WRITE(6,207) ZCOR(1) 

WRITE (8,800) 

FORMAT ('0',* STAGNATION FUNCTIONS') 

WRITE (8,801) 

FORMAT ('0','V,DV,G//D1,D2,D3//BL, BLE, ¥ (1) //DUUDY, DIDY 
*x,CTSR!') 

WRITE (8,207) V,DV,G 

WRITE(8,207) D1,D2,D3 

WRITE (8,207) BL,BLE,2ZA(1) 

WRITE (8,207) DUUDY,DIDY,CTSR 

FORMAT (6F2G. 6) 

pO 4. J=2,8 

DEL=H (J) /™ 

THETA=ABS (ATAN (({ZA (J) -ZA (J-1)) / (R (J) ~R (J-1)))) 
SIT=SIN (THETA) 

CT=COS (THETA) 

JIJ=J 

WRITE (8,205) 

FORMAT ( 8 44% 2 3 28 4 2 3 OR HO TO AIA HH FOI A A IIE HA EE 8 
WRITE (8,202) JJ,X,S (JJ) 

FORMAT ('O',* NEW ELEMENT! ,15,2F20.6) 

DO 1 I1=1,8 

X=X+DEL 

CONTINUE 

CALL FUN (PR,ST,&1,G6,FG,GGG,GK,GL,D1,D2,D3) 
WRITE (8,208) X 

FORMAT (*0', 40K," #k#RRRRERe = X=',F20.6) 
CALL SPFUN(H,U,S,N,X,J3,V,DV) 

RR=R (JJ-1) #DEL*I* (R (JJ) —B (JJ-1)) / (S (JJ) -S (JI 1) ) 
SRR=RR 

IF (NAJEQ.0) RR=1.0 

ZZ=ZA (JJ-1) + 1* (ZA (JJ) -ZA (JI-1)) 78 
DZ2=RR¥RR*FG/V 
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DW 2=RR*RR*¥GGGC/V 
ZS=ZS+ (DZ1+DZ2) *DEL/2.6 
WS=WS+ (DW1+DW2) *DEL/2.9 
Z=ZS/ (RR¥RR) 
W=WS/ (RR¥RR) 
GK 1=Z*DV 
GLi=W*DV 
IF(I.NE.M) GO TO 50 
WRITE (8,204) 
204 FORMAT('O','GK1,Z,Z2S//GL1,W,WS//RM1") 

WRITE(8,207) GK1,2,2Z5 
WRITE (8,207) GL1,W,WS 
IF (2.LT.0.0.08.¥W.L7T.0.6) GO TO 301 
RM1=SQRT (D2* D2*W/ (D3*D3*2Z) ) 
WRITE(8,207) RM1 

50 CONTINUE 
IF (JJ.EQ.2) GO TO 3 
CALL ROOT (PR,ST,R1,G6,GK1, IFLAG,JJ) 
IF(IFLAG.EQ.2) GO TO 5900 

3 CONTINUE 
CALL FUN (PR,ST,R1,G,FG,G6GG,GK,GL,D1,D2,D3) 
DZ1=RR* RR*¥FC/V 
DW 1=RR*¥RR¥*GGG/SV 

EVALUATE PARAMETERS 
THK2=SQRT (2) 
THK3=SQRT (W) 
BL=THK2/D2 
BLE=THK3/D3 
VW=6,9*P/BLE 
vVwi=6.0*P/(R1*BL) 
DUUDY=(2.0+G/6.0) /(BL+BL*PP) 
DTDY=2. 0/ (BLE+BLE* PRP) 
CTSR=DUUDY*V 
IF (I.NE.M) GC TO 51 
WRITE (8,802) 
802 FORMAT('O','V,DV,G//D1,D2,D3//BL,BLE//DUUDY,DIDY,CTSR') 

WRITE (8,207) V,DV,G 
WRITE({8,207) D1,D2,D3 
WRITE (8,207) BL,BLE 
WRITE(8,207) DUUDY,DTDY,CTSR 
WRITE(8,212) VW,VW1 

51 CONTINUE 
DWZ=VW*XCT* DT + NSS¥*CTSR*¥CT*DT 
DWR=VW*SIT*¥DT+NSS*CTSR¥SIT*DT 
DWZN=VW*DT/CT+NSS¥*CTSR*DT/CT 
ZN=ZZ+DWZ 
RN=SRR+DWR 
ZNN=ZZ+DWZ2N 
IF (I.NE.M) GC TO 52 
WRITE (8,209) VW,RN,ZN,ZNN 
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FORMAT ('0','VN,RN,ZN,ZNN!’ ,4F20.6) 
CONTINUE 

IF (I.EQ.M) ZCCR(J)=ZNN 
CONTINUE 

CONTINUE 

RETURN 

IER=1 

WRITE(6,206) IER 

FORMAT ('O*,* ERROR',I4) 
WRITE (6,207) X,G,R1,V,DV 
WRITE (6,207) GK,GL,D1,D2,D3 
WRITE (6,207) 20,WO 

RETURN 

IER=2 

WRITE(6,206) IER 

WRITE (6,207) X,G,R1,V,DV 
WRITE (6,207) GK,GL,D1,D2,D3 
WRITE (6,207) Z2,W,GK1,GL1 
RETURN 

END 


SUBROUTINE BCS (F,X,N,M,GG,DG,H,IER) 
IMPLICIT REAL*4 ({A-H, 0-2) 

DIMENSION F(1),X(1),H (1) ,GG (1) ,DG (1) 
WRITE (6, 100) 

ITER=0 

DO 1 J=2,N 

H (J) =X (J) -X (J-1) 

CONTINUE 

pO 3 J=1,N 

B=2.0 

IF (J.EQ.1) GO TO 4 

IF (J.EQ.N) GC TO 4 
C=H(J3+1) /(H(J) +H (J+1) ) 

A=1.0-C 

D=6. 0* ( (F (J+1) -F (J) ) /H (+1) - (F (5) ~F (J-1)) ZH (9)) 


*/ (H (J) +H (J+1)) 


GO TO 5 
A=0.0 
c=0.0 
D=0.0 
CONTINUE 
Q (J) =H (N+J) U (J) =H (2*N+J) 
P=A*H (N+J-1) +B 
IF (ABS(P).1T.1.E-40) GO TO 10 
H (N+) =-C/P 
H (2#N4+d) = (D-A¥*H (2*N+J-1) ) /P 
GO TO 3 
LTER=130 
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WRITE (6,301) IER 

GO TO 18 

CONTINUE 

DO 141 K=2,)8 

J=N+1-K 

H (2¥*N+J) =H (N+J) *H (2*N+3+1) +H (2*N+J) 
DO 12 I=1,8 

IF (GG (I) .LT. X(1)) GO TO 13 
IF (GG (I) .GT.X(N)) GO TO 14 
DO 15 J=2,N 

JJ=J 

IF (GG (I) .LE.X(J)) GO TO 17 
CONTINUE 

JJ=2 

ITER=131 

WRITE (6,302) IER,I 

GO TO 17 


WRITE (6,302) IER,I 

CONTINUE 

A=X (JJ) -GG (I) 

B=GG (I) -X (JJ-1) 

C=H(JJ) 

GG (I) = (H (2*N+JJ-1) *A¥A*A+H (2*N+JJ) *B¥B*B 
*+ (6. O¥F (JJ-1)-H(2*N+J33-1) *C*C) *A 
2+ (6. 0¥*F (JJ) —H (2*N+Jd) *C¥C) *B) /(6.0*C) 

DG (I) =-1.0*H (2*N+JJ-1) *A¥A/ (2. 0*C) 
*+H (2*N4+3J) ¥B¥B/(2.0*C) + (F (JJ) -F(JJ-1))/C 
% = (H (2*N+JJ) —H (2*N+53~-1) ) *C/6.0 

CONTINUE 

WRITE (6,101) (GG (ZI) ,I=1,M) 

WRITE (6,101) (DG(I) ,I=1,¥) 

CONTINUE 

PORMAT(*0",* OUTPUT FROM CUBIC SPLINE INTERPOLATION’) 
FORMAT (10F9.5) 

FORMAT ("0',//,T10, "ERROR! ,13,//,T10,'NO SOLUTION’) 
FORMAT (*0',//,T10,*ERROR',13,'AT I=",13,//,T10, 
*'SOLUTICN TRIED’) 

RETURN ; 

END 


SUBROUTINE SPFUN(H,F,S,N,X%,J3,SPF,DSF) 
DIMENSION H(1),S(1),F(1) 

A=S (JJ) -X 

B=X-S (JJ-1) 

C=H{(JJ) 

NJJ=N*¥2+5d 

NJJ1=NJJ-1 
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SPF= (H(NJJ1) *A*¥A¥A+H (NDJ) *B¥B*B+t (6. 0¥*F (JJ-1) - 
*H {NJJ1) *C¥C) *A+(6.0*F (JJ) -H (NJJ) *C¥C) *B) / (6.0 ¥*C) 
DSF=-1.0¥*H (NJJ1) *A¥*A/ (2. O*C) +H (NIJ) *B¥B/ (2.0*C) 
*4+ (F (JJ) -F (JJ-1)) /C- (H (NJ J) -H (NJJ 1) ) ¥C/6.0 

RETURN 
END 


SUBROUTINE FUN(PR,ST,R1,G,FG,GG6G,GK,GL,D1,D2, D3) 
PRP=SORT (0. 25+ST/3.0) -0.5 

P=PRP/PR 

PP=P/R1 

R2=R1*R1 

R3=R2*R1 

R4=R3*R1 

D2=(7.4-G/15. 0-G¥G/144 .0+G*PP/10.0+15.6*PP+7, 2*PP*¥PP) 
*/(63.0* (1.04 FF) *¥*2) 

GK=D2*D2*G 

FG=2. O¥D2*((2.0+G/6.0) /(1.04+PP) +6.0*PP) -2.0*G K* 
* (2.0+D1/D2) 

B1=R1-R3/2.0+R4/5.0 
B2=R1/12.0-R2/6.0+R3/8.0-R4/30.0 

B3=R2*2.0-2. 0* R343 2O*R4/5 20 
A1=13.0*R1/15.0-67.0*R3/ 140.047. O*R4/36.0 

A2=0. 8*R1-13.0¥*R3/28. O44. 0*R4/21.9 
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A3=13.0*R1/30.0-13.0*R2/14.,04201.0*R3/280 .0-7. 0*R4/36.4 


A4=0.4*R1-31.0*R82/35.0+39.0*R3/56.0-4. 0*R4/21.8 
A5=13. 0O*R2/7. 0-67. O¥R3/35.0+7.0*R4/12.0 

A6=6 2. O¥ R2/35.0-13.0*R3/7.044, O¥ RU/7.0 

UI= (B1+G*B2+PP¥*B3) /(1.0+PP) 


UTI= (A1+PRP*¥A2Z+G*A3/6.0+G¥*PRP*A4/6,0+PP¥A5+PP¥*PRP* AO) / 


*((1.0+PP)*(1.0+PRP) ) 

D3=U1-UTI 

GL=D3*D3*R2*G 

GGG=2. 0*D3*(2.0/ (PR+PR*¥PRP) +6.0*P) -2.0¥*GL 
RETURN 

END 


SUBROUTINE RCOT(PR,ST,R1,6,GK1,IFLAG, JJ) 
X1=12.0 


G=X1 

CALL FUN (PR, ST,R1,G,FG,GGG,GK,GL,D1,D2,D3) 
FUNC1=GK-GK1 

G=X2 

CALL FUN (PR, ST,R1,G, FG, GGG, GK,GL,D1,D2, D3) 
FUNC2=GK-GK1 
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IF (FUNCI1*FUNC2) 1,2,3 

FUNC1=FUNC2 

X1=X2 

X¥2=X2-4,0 

IF (X2.LT.-12.0) GO TO 490 

GO TO 5 

WRITE (8,101) X1,X2,FUNC1, FUNC2 

FORMAT ('0',* ROOT IS BETWEEN! ,F10.6,*AND',F10.6, 
**FUNCTION VALUES', 2F10.6) 
Xi=(FUNC2*X1-FUNC1*X2) / (FUNC2-FUNC1) 

G=XI 

J=J5+1 

IF (ABS (XI-X1).LT.1.0E-4) GO TO 7 

IF (ABS (XI-X2).1T.1.0E-4) GO TO 7 

IF (ABS (FUNC1-FUNC2) .GT.1.0E-20) GO TO 8 
WRITE (8,102) X1,X2,FUNC1, FUNC2 

FORMAT ('O*,* ROOT IS*,2F10.6,'FUNCTION VALUE*,2F10.6) 
WRITE (8,104) J 

FORMAT ('0',* NUMBER OF ITERATIONS=',16) 
RETURN 

CALL FUN(PR,ST,R1,6,FG,GGG, GK,GL, D1, D2, D3) 
FUNC I=GK-GK1 

IF (FUNCI*FUNC1) 11,12,13 

FUNC1=FUNCI. 

X¥1=xI 

GO TO 6 

FUNC2=FUNCI 

X2=XI 

GO TO 6 

G=X1 

IF (FUNC2.EQ.0.0) G=X2 

IFLAG=90 

WRITE (8, 103) G,X1,X2,I1 FLAG, FUNC1,FUNC2 
FORMAT ('O*,*ROOT IS',3F10.6,' IFLAG=', 14, 
*t FUNCTION *',2F10.6) 

WRITE (8,104) J 

RETURN 

G=XI 

IF (FUNC1.EQ.0.0) G=X1 

IFLAG=1 

WRITE (8, 103) G,X1,X2,1FLAG, FUNC, FUNC2 
WRITE(8, 104) J 

RETURN 

WRITE (6,100) JJ 

FORMAT ('O*',*NC ROOT BETWEEN -12.0 AND 12.0',14) 
IF (GK1.GT.0.0) G=12.0 

IF (GK1.LT.0.0) G=-760 

WRITE (6,105) G 

FORMAT('O*, * ASSUME ROOT=!,F10.2) 
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RETURN 
50 IFLAG=2 
WRITE (6,106) 
106 FORMAT('O',* NUMBER OF ITERATIONS >30') 

STOP 
END 
SUBROUTINE MOVE(ZCOR, NBC, NBCN,X,Y,NN, NDE, NDEN) 
REAL ZCOR(NBCN) ,X (NN) ,¥ (NN) 
INTEGER NBC(NBCN) , NDE(NDEN) 
DIST=ZCOR (1) -4.0 
DO 1 I=1,NBCN,2 
ZCOR (I) =ZCOR (I) -DIST 
NI=NBC (I) 

4 X(NI)=ZC0R(1) 
DO 2 J=2,NDEN,2 
NJ=NDE (J) 
IF (¥ (NJ) .LT.1.15) GO TO 2 
X (NJ) =X (NJ) -DIST 

2 CONTINUE 
RETURN 
END 
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